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Chapter  1 


Introduction 


Flood  flows  in  excess  of  a  reservoir’s  capacity  must  be  passed  down¬ 
stream  in  a  manner  which  does  not  endanger  the  dam  or  surrounding  hydraulic 
structures.  This  is  not  a  trivial  task  as  the  flow  must  fall  a  great  distance  to 
reach  the  river  bed.  These  high  current  velocities  coupled  with  a  free-surface 
can  easily  lead  to  regions  of  low  pressure  in  which  cavitation  may  occur  or  to 
the  formation  of  standing  waves  and  an  uneven  flow  distribution.  Poor  flow 
distribution  will  yield  circulation  and  high  velocities  at  the  base  of  the  spill¬ 
way  (or  outlet  channel)  known  as  the  “stilling  basin,”  resulting  in  downstream 
scour,  potentially  undermining  the  structure  and  causing  both  bank  erosion 
and  stilling  basin  damage. 

Presently,  the  design  of  these  structures  is  accomplished  primarily 
through  large  scale  hydraulic  models.  An  example  is  the  Grapevine  Emer¬ 
gency  Spillway  Model  at  the  Waterways  Experiment  Station  (WES)  shown  in 
Figure  1.1.  A  large  scale  is  necessary  to  eliminate  scale-effects  due  to  the  vis¬ 
cosity  of  the  model  fluid;  i.e.,  the  same  fluid  (water)  is  used  as  in  the  field 
problem.  The  model  shown  is  a  1:60  undistorted  scale  reproduction  of  the 
spillway  section  and  the  apron;  it  includes  a  portion  of  the  upstream  reservoir 
and  the  channel  downstream.  In  the  figure  we  see  the  jump  at  the  toe  of  the 
spillway  apron,  in  which  the  flow  rapidly  decelerates  and  becomes  subcriti- 
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Figure  1.1:  Hydraulic  Spillway  Model. 


cal.  Along  the  apron  the  flow  is  supercritical  (flow  velocity  is  greater  than  the 
wave  speed)  and  in  this  case  is  relatively  smooth.  If  the  upstream  approach 
conditions  are  poor,  oblique  jumps  will  develop  producing  rirculation  in  the 
downstream  subcritical  channel.  This  will  cause  excessive  scour  and  damage. 
The  model  is  used  to  eliminate  the  cause  of  poor  circulation.  The  mode!  is 
rebuilt  or  modified  for  each  design  trial;  e.g.,  by  sidewall  boundary  reshaping 
and  changing  the  bed’s  lateral  curvature.  Each  modification  is  expensive  to 
construct  and  time  consuming.  A  detailed  numerical  model  could  significantly 
reduce  the  design  costs  and  enhance  understanding  of  the  flow  phenomena. 


1.1  Background 

.A  principal  aspect  of  high  velocity  flow  associated  with  spillways  and 
outlet  channels  is  the  formation  of  standing  waves.  A  disturbance  caused  by  an 
obstacle  in  the  flow  or  by  lateral  transition  will  propagate  out  away  from  the 
source.  If  the  flow  velocity  is  faster  than  the  wave  celerity,  c.  the  disturbance 
cannot  travel  upstream.  Instead  it  will  be  swept  downstream  forming  a  wake 
or  standing  wave. 

The  wave  celerity  is  dependent  upon  its  wavelength  [40]  and  is  given 
by 

c  =;  [(5A/2ir)tanh  {27r/i/A)]’  (1.1) 


where,  r  is  the  wave  celerity,  g  is  the  acceleration  of  gravity  A  is  the 
wavelength,  and  h  is  depth.  The  dispersion  of  these  wavelengths  is  a  result  of 
the  reduction  of  the  pressure  gradient  due  to  the  vertical  accelerations  (w-ith 
the  shorter  wavelengths  more  greatly  affected). 
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From  ( 1.1 ),  the  maximum  celerity  is  associated  with  the  longest  wave¬ 
length/depth  ratio  and  has  a  value  Cmox  =  (gh)^  ■  Thus,  if  the  Froude  Number, 
Fr  =  (where,  v  is  the  flow  velocity),  is  greater  than  unity  all  wavelengths 
will  propagate  downstream.  This  is  known  as  supercritical  flow.  Conversely,  if 
F,.  <  1  the  flow  is  subcritical.  A  common  simplification  in  numerical  modeling 
of  these  flows  is  the  long-wave  or  shallow  water  approximation.  This  results 
in  a  hydrostatic  pressure  distribution  with  all  wavelengths  having  the  same 
celerity, 

1.2  Previous  Research 

One  approach  to  flow  over  curved  beds  in  hydraulic  modeling  is  through 
the  use  of  potential  flow  theory.  The  water  surface  is  determined  via  the 
Bernoulli  equation.  An  inverse  approach  in  which  the  coordinates  are  the 
dependent  variables  and  the  stream  function  and  velocity  potential  are  inde¬ 
pendent  variables  has  been  used  by  several  researchers,  e.g.  see  Watters  and 
Street  (62)  for  an  early  application  and  Cassidy  [10]  and  Moayeri  [44]  for  meth¬ 
ods  which  handled  more  general  geometry. 

Examples  of  the  direct  approach  are  given  in  Bettes  and  Bettes  [7] 
and  Ikegawa  and  Washizu  [30].  Here  the  stream  function  and  velocity  potential 
are  the  dependent  variables.  The  grids  adapt  to  follow  the  water  surface.  These 
examples  are  two-dimensional  (longitudinal  and  vertical  resolution)  and  thus 
cannot  determine  the  standing  wave  patterns  and  nonuniform  flow  distributions 
laterally. 

Another  approach  is  to  use  the  shallow  water  equations.  These  equa- 


tions  may  be  derived  following  the  procedure  originated  by  Friedrichs  [24]  and 
extended  by  Keller  [35]  {see  also  Stoker  [56])  by  utilizing  an  asymptotic  expan¬ 
sion  in  a  shallowness  parameter 


■(;y 


where,  I  is  the  radius  of  curvature  of  the  free-surface  and  d  is  a  characteristic 
depth.  To  lowest  order,  the  standard  shallow  water  equations  are  hydrostatic 
and  assume  a  small  channel  slope,  these  equations  are  referred  to  as  the  St. 
Venant  equations  [50].  There  have  been  applications  to  conditions  of  high  ve¬ 
locity  (see  e.g.  [39.  16,  32]),  but  the  assumption  of  mild  slope  and  hydrostatic 
pressure  limits  its  use  in  some  practical  applications  such  as  spillway  simula¬ 
tions. 


A  nonhydrostatic  pressure  distribution  may  be  incorporated  by  in¬ 
cluding  higher-order  terms  (see  Abbott  and  Rodenhuis  [2])  known  as  the  Boussi- 
nesq  terms.  Applications  of  these  equations  to  supercritical  flow  have  been 
made  by  Gharangik  and  Chaudry  [26],  though  the  bed  is  again  assumed  to 
have  a  mild  slope  and  the  equation  set  involves  higher  derivatives. 

Dressier  [20]  produced  a  more  general  set  of  one-dimensional  shallow 
water  equations,  in  which  channel  bed  curvature  is  included  without  resort¬ 
ing  to  incorporation  of  higher-order  terms  in  the  shallowness  expansion.  This 
formulation  leads  to  the  equations 


ht  -f- 


J{hy 


log(J(/t)}  .  ,  log{J(ft)}  1  .  _ 

J{h)K  «  U(h)2^  J{k)  j“  ° 


(1.4) 
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Figure  1.2:  One-dimensional  equation  coordinate  system. 

^  =  0  (1.5) 

J{S3)k  lJ(S3)  J 

P(Si,  S3,  t)  =  pg  COS  e{h-  S3}  +  ^  { J(  A)-*  -  J  (33)"'}  (1.6) 

where  alphabetic  subscripts  indicate  differentiation  and  the  pertinent  variables 
are  (see  Figure  1.2): 

t  is  time 

Si  is  the  coordinate  parallel  to  flow  bed 
S3  is  the  coordinate  orthogonal  to  si 
K  is  the  bed  curvature;  k(si) 

0  is  the  angle  from  horizontal  to  the  tangent  of  the  channel  bed; 

<?(Si) 

h  is  the  depth;  /i(si,<) 

J  is  the  jacobian,  J  =  (\  —  k(si)s3);  J(si,S3) 

P  is  the  pressure 


u  is  the  current  velocity  in  the  si  direction  u  =  ti/J(s3);  u{si,S3,t) 
u  is  the  current  velocity  in  the  direction  at  the  channel  bottom, 
i.  e.,  u  =  u(si,0,t) 

w  is  the  current  velocity  in  the  S3  direction;  u!(sj,S3,<) 
p  is  the  density,  assumed  to  be  a  constant 
ff  is  the  acceleration  of  gravity. 

Dressler’s  expansion  is  applied  to  the  two-dimensional  Euler  equations 
in  the  orthogonal  curvilinear  system  defined  by  Si,S3.  In  the  same  manner  as 
Friedrichs,  irrotationality  is  assumed.  The  other  basic  assumptions  are  constant 
density,  no  lateral  variation  and  that  the  ambient  surface  pressure  is  zero.  The 
usual  kinematic  surface  conditions,  of  no  penetration  at  the  channel  bed  and 
that  a  particle  on  the  free-surface  remains  on  the  surface,  are  enforced.  The 
irrotationality  assumption  is  reasonable  for  converging  flow  [38],  as  is  the  case 
in  the  vicinity  of  the  spillway  crest.  Even  after  the  development  of  the  turbulent 
boundary  layer  this  assumption  for  the  flow  profile  is  quite  reasonable.  The 
resistive  action  of  the  channel  and  eddy  viscosities  resulting  from  turbulence 
can  be  included  in  an  empirical  term  like  the  Chezy  formula.  Furthermore, 
Dressier  [21]  developed  corrections  to  the  Chezy  and  Manning  coefficients  to 
accommodate  bed  curvature. 

Dressier 's  formulation  yields  u,  h,  and  P  to  order  and  w  to  order 
£*.  (The  value  of  w  to  order  is  simply  u;  =  0  which  is  identical  to  the  re¬ 
sult  in  standard  shallow  water  theory.)  Given  the  zero-order  perturbation  for 
the  solution  u  and  h,  the  next  higher  approximation  of  w  may  be  calculated. 
This  treatment  utilizes  evaluation  of  the  equations  at  the  channel  bed  to  re- 
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move  some  complicating  terms.  Sivakumaran  (51,  52j  generalizes  the  derivation 
further  to  a  two-dimensional  surface.  Again,  irrotationality  is  assumed.  The 
one-dimensional  equations  when  evaluated  at  the  bed  then  become: 


J{h)ht  +  ^ 

US\ 


~-\og{J(h)] 


=  0 


+  gE,,  =  0 


(1.7) 

(1.8) 


where, 


E{si,t)  2  C-f-Acos^-f  — -F  ^ 


•2 

tc* 


P-Po 

PQ 


pg  J{hy2g 

=  [A  -  sajcosd-t- [j(A)~*  -  J(s3) 
I  d  *■ 


-2 


{P 


2g 


tv  = 


J (^3) 


u 


log{J(53)} 


(1-9) 


(1.10) 


and  C  is  the  elevation  of  the  bed  above  some  reference  level  with  Po  the  ambient 
pressure  at  the  free-surface.  Both  of  these  sets  shall  be  referred  to  as  Dressier ’s 
equations  throughout  the  present  study  as  these  two  formulations  are  equiva¬ 
lent.  This  model  appears  applicable  for  —0.85  <  nh  <  0.5  and  Sivakumaran 
(51,  53]  has  demonstrated  a  good  agreement  with  experimental  results  over  the 
wider  range  — 2  <  «h  <  0.54. 


1.3  Objectives 

The  principal  objective  of  this  present  research  is  to  develop  and  eval¬ 
uate  a  generalized  set  of  shallow  water  equations  containing  bed  curvature 
effects  to  simulate  important  aspects  of  flow  over  curved  beds,  specifically  as¬ 
sociated  with  spillways  or  outlet  works. 

Subordinate  objectives  in  support  of  the  primary  aim  are: 
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1.  Via  a  perturbation  analysis  develop  the  generalized  set  of  shallow  equa¬ 
tions  without  excluding  vorticity  about  the  direction  normal  to  the  bed. 

2.  Develop  a  one-dimensional  finite  element  model  for  preliminary  tests  of 
these  equations  with  a  particular  emphasis  on  the  effects  of  curvature. 
Results  will  be  compared  to  the  steep-slope  shallow  water  equations  which 
do  not  include  bed  curvature. 

3.  Develop  a  two-dimensional  numerical  model  of  this  equation  set. 

(a)  Make  self-consistency  test  of  the  model  to  demonstrate  that  the 
model  matches  the  derived  equations. 

(b)  Test  the  numerical  method  employed  by  comparison  to  flume  water 
surface  data. 

(c)  Make  comparisons  to  flume  results  for  these  equations,  the  standard 
steep-slope  shallow  water  equations  and  the  St.  Venant  equations. 

1.4  Important  Results 

The  important  results  from  this  work  are: 

1 .  The  development  of  a  system  of  nonhydrostatic  two-dimensional  shallow 
water  equations  that  are  not  restricted  to  irrotationality. 

2.  A  comparison  with  the  standard  steep-slope  shallow  water  equations  re¬ 
veals  unexpected  errors  upstream  of  a  simulated  dam  crest, 

3.  A  numerical  scheme  that  shows  promise  for  super-  and  subcritical  flow 
in  hydrodynamic  modeling. 
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4.  An  in-depth  comparison  of  the  usefulness  and  limitations  of  these  new 
equations  along  with  more  traditional  St.  Venant  and  steep-slope  stan¬ 
dard  shallow  water  equations. 

5.  A  numerical  model  capable  of  numerically  simulating  hydraulic  flow  along 
chute  spillways  and  outlet  works. 

6.  A  data  set  of  water  surface  elevations  from  an  outlet  works  that  is  of 
general  use. 

1.5  Outline  of  Treatment 

A  general  set  of  shallow  water  equations  which  is  nonhydrostatic  and 
which  allows  vorticity  is  derived  in  Chapter  2.  The  derivation  involves  an 
asymptotic  expansion  in  a  shallowness  parameter  for  the  bed-fitted  Euler  equa¬ 
tions.  Preliminary  testing  of  these  equations  in  one  dimension  is  conducted  in 
Chapter  3.  A  comparison  with  a  standard  steep-slope  shallow  water  equation 
set  is  made.  In  Chapter  4  we  detail  the  development  of  the  two-dimensional 
model  with  special  attention  to  a  numerical  scheme  to  better  treat  the  highly 
convective  flow  conditions.  A  more  in-depth  set  of  tests  is  then  conducted  for 
the  two-dimensional  generalized  shallow  water  equations  (Chapter  5).  We  first 
test  the  validity  of  the  code  for  several  test  problems.  The  numerical  scheme 
and  the  treatment  of  boundary  conditions  are  then  compared  with  the  mea¬ 
sured  water  surface  in  a  flume  for  a  supercritical  transition.  The  influence  of 
curvature  and  the  mild-slope  assumption  (for  the  St.  Venant  equations)  are 
made  apparent  by  comparison  with  previously  collected  flume  data.  The  most 
general  test  of  the  equations  is  a  comparison  to  water  surface  of  an  outlet- 
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works  fiume.  Here  the  bed  contains  curvature  and  there  are  lateral  transitions. 
Tests  on  this  configuration  are  made  for  the  St.  Venant  equations,  the  steep- 
slope  standard  shallow  water  equations  and  the  generalized  set  by  comparison 
of  predicted  water  surface  elevations  with  flume  results.  Chapter  6  contains 
conclusions  concerning  the  usefulness  of  these  equations,  their  advantages  and 
disadvantages,  and  recommendations  for  additional  research. 


Chapter  2 


Equation  Development 


The  free-surface  and  the  nonhydrostatic  pressure  aspects  of  the  flow 
were  included  in  the  previous  equations,  but  vorticity  about  the  bed- normal  di¬ 
rection  was  excluded.  It  is  important  that  eddy  patterns  and  vorticity  resulting 
from  sidewalls  and  the  drag  due  to  bottom  friction  be  reflected  in  the  equations 
to  distinguish  the  design  alternatives.  Also,  the  variation  in  pressure  due  to 
curvature,  in  fact,  may  rival  that  of  the  hydrostatic  pressure.  These  problems 
cannot  be  adequately  modeled  using  the  previous  perturbation  formulations. 
The  present  approach  yields  a  more  general  formulation  that  does  not  restrict 
vorticity  about  bed-normal  axes  while  including  bed-curvature  effects. 

The  derivation  developed  here  employs  concepts  common  to  the  stud¬ 
ies  of  Friedrichs,  Keller  and  Dressier  and  also  involves  easing  of  the  irrotational- 
ity  restriction  with  extension  to  a  two-dimensional  surface.  The  basic  approach 
is  to  use  an  asymptotic  expansion  of  the  dependent  variables  of  the  three- 
dimensional  Euler  equations  (written  for  an  orthogonal  curvilinear  coordinate 
system)  in  a  shallowness  parameter  e. 

2.1  Basic  Equations 

The  derivation  begins  with  the  Euler  equations; 

V-v  =  0  (2.1) 
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Vt  +  (v  •  v)  —  V  X  u?  +  -VP  —  g  =  0  (2.2) 

2  p 

where,  u;  is  the  vorticity  vector,  and  g  is  the  body  force  per  unit  mass.  The  free- 
surface  kinematic  boundary  condition  requires  that  a  particle  on  the  surface 
remain  on  the  surface,  so  that 

/i(si,S2,t)  -  S3  =  0  (2.3) 

This  may  be  written 

hi  +  v(si,S2,h,tyVh  —  tv{si,S2,h,t)  (2.4) 

where,  u,  v,  and  w  are  in  the  sj,  S2,  and  S3  directions,  respectively.  There  is 
also  a  bottom  kinematic  boundary  condition  that  the  velocity  normal  to  the 
bed  be  zero.  This  implies 

U)(5i,S2,0,f)  =  0  (2.5) 

The  pressure  at  the  free-surface  is  a  constant,  here  taken  as  the  reference  value 
zero, 

P(si,S2,/i,t)  =  0  (2.6) 

The  irrotationality  condition  implies  that  the  vorticity  vector  components  in 
the  Si  and  S2  directions  are  zero. 


—  0  (2.7) 

The  coordinate  system  shown  in  Figure  2.1  is  a  mutually  orthogonal 
system.  The  coordinate  directions  si  and  S2  are  in  fact  curvilinear,  and  S3 
is  normal  to  the  surface  so  defined.  Therefore,  53(1,  y,z)  =  c  (where  i,y,2 
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Figure  2.1:  Bed-fitted  coordinate  system. 

are  cartesian  coordinates)  defines  a  coordinate  surface  above  the  bed.  An 
infinitesimal  vector  dg  that  lies  within  the  bed  has  a  length  given  by 


where, 


Ci  = 

1  dx 

Wi' 

1  dx 

dx\^ 

C2  — 

032) 

and  X  is  a  position  vector.  For  the  surface- normal  coordinate  =  1- 

As  one  moves  along  the  S3  direction,  (i  and  Cj  are  scaled  in  relation  to  S3 
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as  ji  =  Ci(l  —  and  similarly  for  the  S2  direction.  This  formulation  of 

the  metrics  requires  that  si  and  sj  not  only  be  orthogonal  but  also  principal 
directions.  The  curvature  in  the  Si  direction,  for  example,  is  given  by,  (see 


[41]): 


Ki  = 


dx 

9ai 


931 


/9x..^\2  /  9x  9x  3 

V^ai  931/  \9»i  ) 


(2.8) 


Where  N  is  the  unit  vector  normal  to  the  s\  —  surface;  if  is  parallel  to 
the  Si  tangent  vector  then  the  curvature  along  si  is  an  extremum  and  si  is  a 
principal  direction.  The  normal  vector  remains  in  the  si  —  S3  plane;  s^  is  normal 
to  Si  and  hence  is  also  a  principal  direction.  Introducing  this  coordinate  system 
(2.1)-{2.7)  and  simplifying,  we  find  the  following  equation  set. 


Continuity  Equation; 


djjiv)  djjijzw) 
ds\  9s 2  dss 


(2.9) 


Si  Momentum  Equation: 


u  V 

u,  -f  r-u,,  +  — +  wu, 
Ji  J2 


^  «t>  djj 
jl]2  9si  Jij2  dsj 


S2  Momentum  Equation; 


u  V 

Vt  +  — +  wv,, 
Ji  h 


J1J2  ds2  J1J2  9si 


vw  +  - —  P,J  -  92  =  0 

J2P 

(2.11) 


S3  Momentum  Equation; 


u  V 

tl’t  +  —  +  —  lUjj  +  WW,J  + 

Jl  J2 


v'  +  -/"33 -53  =  0  (2.12) 
P 
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Irrotationality  Condition; 


W,,  -  (;2V)„  =  0 

(2.13) 

=  0 

(2.14) 

Free-Surface  Boundary  Conditions  (at  S3  =  h  ): 

ht  +  =  u>  (2.15) 

Jl  J2 

P{s^,S2,h,i)  =  0  (2.16) 

Channel  Bed  Boundary  Condition: 

jz;(.s,,jr2,0,0  =  0  (2.17) 

These  equations  are  quite  complex  and  our  objective  now  is  to  extract  the 
important  simplifications  of  these  equations  that  result  when  the  flow  depth  is 
small. 

In  the  manner  of  Friedrichs,  the  equations  may  be  nondimensionalized 
and  the  dependent  variables  expanded  in  powers  of  e.  Typical  length  scales 
are  the  depth  and  the  free-surface  radius  of  curvature  As  in  shallow 
water  theory  the  relationship  with  velocity  and  thus  time  is  assumed  by  the 
approximate  celerity  of  a  free-surface  wave,  (gd)^  [56]. 

Nondimensionalizing,  we  let  a,  7  and  r  be  the  new  independent 
variables  so  that 

s,  = /q  u  =  igd)^  ii  (1  =  Cl  K,  ==  Ki/d  Jl  =Ct(l-^i/?) 

$2  =  l-y  V  =  (gd)^  V  (2  =  C2  «2  =  •^2/d  J2  =(2(1-  ^2/^) 

S3  =  dd  XV  =  P  =  pgdw  t  —  — ‘-jT  =  £ 

h  =dY 


(2.18) 
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where  the  tildes  ( ‘  )  indicate  nondimensional  quantities.  Recasting  the  gov¬ 
erning  equations  in  the  dimensionless  form  and  then  dropping  the  tildes  for 
convenience,  equations  (2.9)-(2.17)  become 


Continuity  Equation: 


£  {(jju)a  +  0>K}  +  -  0 


(2.19) 


a  Momentum  Equation: 


{  u  V  djj  uv  dj,  r. 

<  Ur  +  — Uo  +  — U.,  -  -r— -5—  +  --r~  +  - - F\  WU0  -  —  UU>  =  0 

I  Ji  J2  JiJiOa  Jij2d'r  Ji  J  V-^/i 


7  Momentum  Equation: 


f  ,  ^  djl  ,  dji  ,  TTy  _ 

<  Ur  +  -rVa  +  —Vy - r-: — r - H  - hr - G  >  -f 

I  Ji  32  3x32  07  3x32  da  J2  J 


(2.20) 


WV0  —  1  —  uu;  =  0 

V  ;  Jj 

(2.21) 


0  Momentum  Equation: 


|u;r  +  +  yWy  +  -f  -  //|  +  WWq  =  0  (2.22) 


Irrotationality  Condition; 


-  {J»0  =  0 


(2.23) 


(jiu)o  -  u;o  =  0 


(2.24) 


Free-Surface  Boundary  Conditions  (  at  /?  =  K  ); 


|yr  +  -y;  +  -n}  = 

I  3i  32  j 


Tr(a,j,y,T)  =  0 


(2.25) 


(2.26) 


Channel  Bed  Boundary  Condition: 


u;(q,7,0,t)  =  0 


(2.27) 
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2.2  Perturbation  Analysis 

A  perturbation  analysis  in  e  is  now  developed  for  the  flow  equations 
in  the  curvilinear  system.  We  expand  each  dependent  variable  in  a  power  series 
in  e  as  follows: 


u(a,7,/3,r;£:)  =  u<°'(a,7,/3,T)  +  eu<'>(a,7,/3,T)  +  --- 
t;(a,7,/3,r;£)  ?=  u<®>(a,7,/3,T)-fev<’\Q,7,;5,T)  +  --- 
ti;(o;,7,/?,T;£:)  =  t/;<°^(a,7,y?,  r)  +  £:i;<*){a,7,/?,  r)  +  •  •  •  (2.28) 

7r(a,7,;3,T;£)  =  7r<°^(Q,7,/3,T)  +  £t<‘>(c»,7,/?,t)  +  •  •  • 

K(a,7,^,T;£)  =  V'<°>(a,7,/?,T)  +  £r<*^(a,7,/3,r)  +  •  •  • 


These  expansions  are  substituted  into  equations;  (2.19)-(2.27)  yielding: 

£{L;2u(a,7,^,T;£)]^  +  [jiv{a,^,0,T;€)\^} 

+  [jU2UJ(a,7,/?,r;£)]^  =0  (2.29) 

£  |[u(a,  7, 0,  r;  e))^  +  j-«{o,  7,  /3,  r;£)  [u(a,  7, 0,  r;  e)]^ 

+— v(o,  7,  /?,  r;  £)  {u(o,  7,  /?,  r;  £)] 

J2 


4-4-l7r(a,7,/3,T;£)] 


,  1  -  ,3 

+  u?(a,  7,  /?,  r;  £ )  (u(q, 7,  t;  £ ))g 


— i_[v(a,7,/3,T;£)]^^  -f  -^u{a,7,/?,T;e)v(Q,7,/3,r;£)~ 
J1J2  aa  ;ij2  07 


-I 


J 


u(a.  7,  /?,  r;  £)u;(a,  7,  /9,  t;  £)  =  0 


1 


(2.30) 


£  |[v(Q',7,/3,r;£)]^  +  j-u(Q,7,i3,T;£)  b(o,7,/3.T;£)]^ 
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32 

dj 


32 
1 

3i32 

+  4-  [7r(a,'7,  /?,  t;  e)]^  —  cl  +  w{a,  7, 0,  t;£)  [v(q:,  t\ f)]^ 
32  ) 


-  u(a,7,/?,r;£)u)(or,7,^,T;£)  =  0 

e|[io(a,7,i3,T;£)]^  +  ju{a,f,0,T\e)[w{a,f,0,T\€)]^ 

+  iu(o,  7, 0,  r;  £)  (uj(a,  7, 0,  r;  e)]^ 

32 

+  [u(Q,7,/3,T;£)f  +  (u(a,7,^,T;£)f 

+  (?r(o:,7,i3,T;£)]^- 

+i(7(a,  7,  /5,  t; £)  Ha, 7, 0,  r\ e))^  =  0 

(u;(a,7,/3,T;£)]^  -  \j2v{an-,0^'^\^)\o  =  0 

[jiTi(Q,7,/?,T;£)l^  -  HQ,7,/?,'r;£)]o  =  0 

£{[y(Q,7,^,r;e)U 

+  -r-; - ^:r-T«(0’7i^^'r;£)l>'(a,7.^>T-;£)]„ 

ji(a,7,K,T) 

+  — ; - ^-77— 7.  C)  7, 0^ 

32{a,-t,Y,T)  ) 

=  u){q,7,V',t;£) 

7r(a,7,K,  t;£)  =  0 
u;(q,7,0,t;£)  =  0 


(2.31) 


(2.32) 

(2.33) 

(2.34) 


(2.35) 

(2.36) 

(2.37) 
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Perturbation  equations  are  developed  by  collecting  terms  of  identical 
order  in  e.  Thus,  these  equations  reflect  the  relative  significance  of  the  flow 
shallowness.  Beginning  with  the  lowest  order  effect  {e°),  the  continuity  equation 
in  conjunction  with  the  channel  bed  boundary  condition  imply  that 

=  0  (2.38) 


The  irrotationality  conditions  yield  the  relationship  between  the  flow  velocities 
and  depth 


„(o)  ^  (4i 

ii 

where, 


(2.39) 

(2.40) 


V  =  (2.41) 

u  =  u^°’(o,7,0.  r)  (2.42) 


Now  consider  the  first-order  perturbation  effect.  We  obtain  the  continuity 
equation  contribution 


Note  that  this  can  then  be  integrated  with  respect  to  0  over  the  depth  to  yield; 

d  ftiO  r  /«2  -  «! 


u;'^>(r)  = 


1 

■  d 

MYwy) 

da 

d  ■ 

^7 

Kl 


«2 


4.  +  yMIlylo) 

/i(K)  "  /2(K)  " 


(2.43) 
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where, 

fi  =  1  - 

/2  (K<°‘)  =  1  -  K2K<"> 

Substituting  in  the  free-surface  equation  resulting  from  the  terms  of  order  e  we 
obtain 


=  0  (2«) 

Now,  the  momentum  equation  for  the  normal  direction  when  integrated  with 
respect  to  0  yields  the  zero-order  perturbation  pressure  contribution, 

^  {i^  ((/.  (K'”'))-’ -  (/.  (/?)r) 

+  ((/,  (K'”'))'"  -  (A  («)-’)  }  +  H  {r*"'  -  /j}(2.45) 

The  terms  within  the  first  set  of  braces  are  the  contribution  to  the  pressure  from 
the  curvature  “centrifugal”  effect  and  the  second  set  is  simply  the  hydrostatic 


component.  The  momentum  equation  in  the  q  direction  can  be  simplified  by 
the  realization  that  =  0  due  to  (2.40).  This  implies  then, 


Ji  32  3x32  da  >1>2  d-f  ji 
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2.3  Generalized  Shallow  Water  Equations 


We  now  transform  back  to  dimensional  form  and  drop  the  pertur¬ 
bation  index  for  convenience  and  find  the  equations  given  below.  The  mass 
conservation  equation  is  derived  from  equation  (2.44),  the  monnentum  equa¬ 
tions  from  equations  (2.46  and  2.47),  the  pressure  from  equation  (2.45),  and 
the  vertical  velocity  from  equation  (2.43). 


Ji{h)j2{h)h,  +  -~ 
d 


.<,2  r 

u—  ( 

/Kj  -Kj\ 

\  Ki  J 

+ 


dS2 


L  *  '  '  * 

d32  ,  uv  dji  ,1  „ 


u  V  dj2  uv  dji  I  ^ 

Ut  +  -T-Uj,  +  - - ^  ~ — ^  5i  ==  0 

Jl  J2  J1J2OS2  JiP 

U  V  dji  uv  032  1 

Ut  +  -h  -  77-^  +  T-7-”  +  -T— -  52  =  0 

J\  h  J1J2OS2  J132OS1  32P 

— — -  K1S3)  -f  K2S3 

> 


w 


>(53)  =  “ 


1 


JIJ2 


_d 

ds 

d 

ds 


:{?(( 

{5(t 


«1 

«1  -  «2 
K2 


[  \1  —  K2/1/  \1  —  K2S3 


-)log{l  - 

log(l  -  K2S3)  +  K1S3 


+  553(^1  — 


ti(53)  =  7-^ 

1  —  ^1^3 

^(33)  =  :j - - - 

1  -  K2S3 

Ji  =  Cl  (1  ~ 


(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 

(2.53) 

(2.54) 

(2.55) 


^2  —  C2  (1  —  «^2'S3) 


(2.56) 
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-  Kth) 

(2.57) 

i2(fe)  =  C2(l- 

-  Kih) 

(2.58) 

This  system  of  equations  can  be  used  to  solve  for  h,  u,  v,  and  P  and,  from 
this  information,  w. 

These  equations  represent  the  zero-order  approximation  for  /i,  u,  and 
u;  and  as  such  are  only  accurate  when  the  flow  is  shallow  relative  to  longitudinal 
and  lateral  features.  This  means  that  the  water  surface  reasonably  parallels 
the  bed.  This  is  true  of  all  shallow  water  equations.  The  bed-normal  velocity 
(equation  2.52)  is  the  first-order  representation  of  u;,  (the  zero-order  is  that 
the  bed  normal  velocity  is  identically  zero).  This  is  the  condition  reflected 
in  the  solutions  of  A,  u,  and  v.  The  first-order  approximation  for  w  results 
from  inferences  drawn  using  the  mass  conservation  equation.  Its  effects  are  not 
“felt”  in  these  shallow  water  equations.  Higher-order  expansion  terms  would  be 
necessary.  It  is  important  to  keep  in  mind  the  limitations  on  these  generalized 
shallow  water  equations. 


Chapter  3 


One-Dimensional  Problem  and  Exploratory  Studies 

3.1  Introduction 

Shallow  water  theory  has  been  extensively  studied  both  analytically 
and  numerically  for  a  wide  class  of  free-surface  flows.  However,  there  have  been 
relatively  fewer  shallow  water  studies  of  free-surface  flows  over  curved  beds  such 
as  those  encountered  in  spillway  designs.  The  curved  bed  case  presents  certain 
diflBculties  both  analytically  and  numerically  when  compared  with  the  flows 
usually  encountered  in  shallow  water  calculations,  since  the  pressure  distribu¬ 
tion  is  now  decidedly  not  hydrostatic.  As  an  example,  in  flow  over  a  spillway, 
the  bed  curvature  produces  an  inertial  acceleration  and  thereby  an  apparent 
force  which  is  comparable  in  size  to  the  hydrostatic  pressure.  The  spillway  ca¬ 
pacity  therefore  is  strongly  influenced  by  these  curvature  effects.  In  this  chapter 
we  first  study  the  properties  of  the  one-dimensional  problem  and  develop  a  cor¬ 
responding  finite  element  treatment.  The  associated  numerical  studies  permit 
a  preliminary  study  of  the  overall  approach.  This  in  turn  is  used  to  guide  the 
development  of  the  final  two-dimensional  model  and  analysis  in  the  following 
chapters. 

The  main  thrust  of  the  present  study  is  directed  towards  analysis  and 
finite  element  approximation  of  this  problem,  including  the  influence  of  the 
curved  bed.  The  analysis  extends  the  approach  of  Dressier  (20]  by  means  of 
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a  perturbation  analysis  leading  to  a  new  formulation  of  the  problem  (recall 
Chapter  2)  and  a  generalized  set  of  shallow  water  equations  which  includes 
bed  curvature  effects.  In  the  limit  of  zero  curvature  the  classical  shallow  water 
equations  are  recovered.  The  problem  is  complicated  by  the  fact  that  there  are 
possible  transitions  between  subcritical  and  supercritical  flows  depending  upon 
the  local  Froude  number  [27].  A  further  issue  arises  in  the  representation  of 
the  hydraulic  jump  which  forms  near  the  base  of  the  spillway.  This  hydraulic 
jump  enters  as  a  discontinuity  in  the  solution  of  the  mathematical  problem. 

As  with  standard  shallow  water  theory,  certain  simplifying  assump¬ 
tions  also  arise  here  implicitly  from  the  perturbation  analysis  and  the  neglect 
of  higher-order  terms.  The  primary  one  is  that  accelerations  normal  to  the 
bed  are  nonexistent  and  thus  the  dispersive  character  of  shorter  wavelengths 
(in  which  wave  speed  depends  upon  wavelength)  will  not  be  represented.  The 
resulting  flow  variables  can  only  comprehend  longer  wavelengths  in  which  wave 
speed  depends  upon  water  depth  and,  in  this  particular  set  of  equations,  bed 
curvature.  The  overall  features  of  the  hydraulic  jump  can  be  captured  as  weak 
solutions  containing  a  discontinuity.  However,  once  again,  vertical  accelera¬ 
tions  and  short  period  waves  are  not  produced.  A  discussion  of  the  equations 
and  details  of  the  jump  representation  are  included  in  Sections  3.2  and  3.3. 

In  the  present  treatment,  an  approximate  formulation  based  on  finite 
elements  (Section  3.4)  is  constructed  for  the  dominant  perturbation  terms  in 
the  solution  to  the  curved  bed  problem.  To  accommodate  the  wave-like  be¬ 
havior  and  produce  a  stable  method,  a  special  artificial  viscosity  formulation  is 
employed.  This  form  is  motivated  by  the  eigenvalue  properties  of  the  system 
and  the  existence  of  subcritical  and  supercritical  regions.  The  resulting  scheme 


26 


introduces  numerical  dissipation  in  the  discrete  model.  Partly  as  a  consequence 
of  this,  the  hydraulic  jump  is  approximated  over  several  elements  but  still  does 
yield  good  agreement.  Numerical  experiments  are  conducted  for  a  represen¬ 
tative  problem  corresponding  to  shallow  water  flow  over  a  “bump”  (Section 
3.5). 

3.2  One-Dimensional  Equations 


The  perturbation  analysis  implies  for  mass  conservation  (by  depth 
integration  of  (2.49)  ) 


^ _ 9  j'plogj(ft)]  ^ 

dt  j{h)  9s  \  Kh  ) 


and  for  conservation  of  momentum  (depth- integrated  (2.50)) 


MM  ^  [j(/l)J 


+  5  gh^  cos  6 


Kuwdsj  -f  j{hf2)ghsm0  -1- 


g»^|p|p 

Cj[h/2)hy^ 


=  0 


(3.1) 


(3.2) 


where  the  Manning  equation  is  utilized  to  account  for  friction  and  turbulence 
losses.  Here,  for  notational  convenience,  we  have  set  j  for  j\  and  s  for  si 
with  p  =  uh  where  «  is  the  velocity  at  the  bed  surface,  S3  =  0.  Further,  d  is 
the  angle  between  the  bed  and  the  horizontal  direction  (see  Figure  1.2).  n  is 
Manning’s  bed  friction  coefficient,  and  C  is  a  conversion  factor  needed  since  n 
is  not  dimensionless,  (1.0  for  metric  units  and  2.21  for  English  units). 


Solution  variables  u  and  w  now  denote  the  zero-order  and  first-order 
perturbation  contributions  for  the  tangential  and  normal  velocity  components, 
respectively  (the  term  w  arises  in  the  depth  integration  and  can  be  written 
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in  terms  of  zero-order  h  and  u,  this  is  simply  a  convenience).  The  integral  m 
equation  (3.2)  can  be  evaluated  directly  by  substituting  the  known  relation¬ 
ship  for  u  over  depth  and  including  the  relationship  (2.52)  for  the  first-order 

contribution  to  w.  We  have 

,  f  log  j(<i)  + /  p  ^  _  1  ^  (3.3) 

Kuwdss  =  I — ;5(Ap- J  \ « y  a,  I  2  [Hk)\  as 


Since 


''®o{ 
lim  ~\ 

*->0  K 


K— K 

log  jjh)  -f  Kh 
Kj{h) 

log  j{h)  -h  Kh 
Kh 


\  - 
}  =  “ 
]  = 


it  follows  that  the  terms  involving  curvature  k  in  the  denominator  of  (3.3)  must 
remain  bounded  in  the  limit  as  k  approaches  zero. 


3.3  Discussion 


Observe  that  the  pressure  force  is  represented  in  equation  (3.2)  by 


•  1 2 

-K  — ^  +  \gh^  cos  0 

2  [j(/i)J 


where  the  first  term  can  be  associated  with  the  centrifugal  effects  and  the 
second  is  due  to  the  hydrostatic  pressure.  Th^c.  the  pressure  distribution  is 
indeed  not  hydrostatic. 

Also  note  that  if  k  =  0,  these  flow  perturbation  equations  (3.1)  and 
(3.2)  reduce  to 

^  +  ^  =  0  (3.7) 

dt  ds 
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^  Jp*  .  I  i.2  A  ^  u  ■  a  ^  ?"^lplp  n 

which  correspond  to  the  shallow  water  equations  in  conservation  form  for  the 
case  where  the  channel  slope  is  significant.  If  channel  slope  is  not  significant 
we  recover  the  classical  shallow  water  equations. 


The  characteristic  differential  equations  for  the  hyperbolic  system 
(3.1  H3. 2)  are 

(3.9) 


ds 

di 


u 


mv 


with 


log  j{h) 


This  implies  that  in  this  model  wavespeed  is  independent  of  wavelength.  That 
is,  all  wavelengths  travel  at  the  same  speed,  which  depends  on  the  depth,  grav¬ 
itational  and  centripetal  acceleration.  A  region  of  flow  is  considered  subcritical 
if  <  and  if  it  is  supercritical.  The  set  of  equations 

describes  flow  in  these  regimes.  However,  the  transition  from  supercritical  to 
subcritical  must  occur  through  a  hydraulic  jump.  This  may  be  mathemati¬ 
cally  characterized  by  a  weak  solution  of  (3.1)-(3.2)  containing  a  discontinuity. 
The  flow  variables  will  not  allow  the  production  of  vertical  acceleration  and 
short  period  waves  associated  with  a  jump,  but  the  overall  flow  features  can  be 
preserved. 


The  equations  developed  here  model  the  total  momentum  flux  through 
a  cross-section  whereas  Dressier  considers  momentum  flux  per  unit  depth.  In 
one  dimension  his  approach  corresponds  to  conserving  mechanical  energy.  How¬ 
ever,  in  reality  steep  gradients  in  water  surface  elevation  arise  in  the  form  of 
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hydraulic  jumps  where  (physically)  mechanical  energy  is  transformed  to  other 
forms  of  energy  through  flow  turbulence,  generation  of  short-waves,  etc.  [Ij; 
this  is  apparent  in  Figure  1.1.  Hence  energy  will  not  be  conserved  but  momen¬ 
tum  is  conserved.  In  the  mathematical  treatment  of  the  hydraulic  jump  the 
downstream  (tailwater)  elevation  is  prescribed  and  the  preceeding  observations 
then  imply  that  the  differences  in  the  present  formulation  to  that  of  Dressier 
would  be  manifest  through  the  location  of  the  hydraulic  jump.  The  new  model 
should  predict  a  hydraulic  jump  location  that  is  closer  to  the  physical  situa¬ 
tion.  Since  the  hydraulic  jump  is  an  important  feature  of  the  flow  from  the 
standpoint  of  the  engineering  design  problem,  this  distinction  is  important. 


Introducing  the  weighted  residual  projection  in  (3.1)  and  (3.2)  against 
test  function  4)  and  integrating  by  parts  with  respect  to  s ,  we  obtain 


f  f 


dt 


ds 


j  ds  — 


9<P  Un  =  0 


and 


x( 


dp  d<t> 

u  \dt'^  ds 


-7^0  -  r—  -f  /(s,<,  /i,p)<A  —  d  Kuwds^  ds  +  r({>  jan  =  0 


where  fl  is  the  domain,  50  is  the  boundary  of  0,  and 

P  logXA) 


T  = 


nh 

-^  +  u 

hj{k)  2  [j{h) 


■f  cos  0 


(3.10) 


(3.11) 


(3.12) 


The  weak  form  of  the  transient  problem  (3.10)  -  (3.11)  admits  discontinuous 
solutions  that  model  the  actual  hydraulic  jump  when  the  shallow  water  flow 
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encounters  the  deeper  water  in  the  “stilling  basin”  downstream.  Mathemati¬ 
cally  this  discontinuity  arises  from  the  “hard”  downstream  condition  where  the 
tailwater  elevation  is  specified.  Thus  the  solution  to  the  weak  problem  may 
have  a  jump  discontinuity  in  p  and  in  h. 

Since  we  are  primarily  interested  in  the  steady  flow  situation,  let  us 
consider  the  steady  form  of  (3.10)  -  (3.11)  for  an  arbitrarily  small  region  R 
containing  discontinuity  point  s  with  test  function  0  zero  on  the  boundary  of 
R  and  outside  R.  Then,  for  conservation  of  mass  (3.10)  implies,  in  the  limit 
as  R  0,  M  =  o  (where  [qj  denotes  the  jump  in  mass  flux  q  at  s).  That 
is,  mass  is  conserved  through  the  hydraulic  jump.  For  the  transient  problem  a 
similar  treatment  applying  the  divergence  theorem  in  (s,t)  yields  the  relation 
=  W  for  the  speed  of  propagation  of  the  discontinuity. 

In  the  case  of  zero  curvature,  k  =  0,  the  governing  perturbation 
equations  reduce  to  the  familar  shallow  water  equations.  A  similar  analysis 
to  that  given  above  can  then  be  applied  to  the  momentum  equation  to  yield 
[r]  =  0  for  the  steady-state  problem  and  \p]^  —  [cj  for  the  transient  problem. 
When  «  ^  0,  the  presence  of  the  term  nuwdsy  in  (3.11)  now  provides  a 
non-zero  contribution  for  (rj  at  s.  In  typical  spillway  applications,  the  bed  is 
curved  but  the  hydraulic  jump  occurs  downstream  in  a  region  of  zero  curvature. 

3.4  Finite  Element  Approximation 

3.4.1  Choice  of  Basis 

Next,  we  consider  approximate  analysis  of  the  weak  statements  in 
(3.10)  and  (3.11)  on  a  discretization  of  finite  elements.  The  underlying  contin- 


31 


uous  problem  is  comprised  of  a  first-order  hyperbolic  system  in  elevation  and 
velocity  variables.  Discretization  may  give  rise  to  spurious  oscillatory  modes 
in  the  shallow  water  problem.  The  presence  and  nature  of  spurious  modes  of 
wavelength  2As,  where  As  is  the  mesh  spacing,  has  been  investigated  previ¬ 
ously  [46,  59).  It  has  been  shown  for  the  linearized  shallow  water  equations 
that  piecewise-linear  approximation  of  u  with  piecewise  constant  h  does  not 
generate  an  oscillatory  mode.  However,  in  the  present  analysis,  the  curvature 
terms  preclude  a  piecewise-constant  basis.  A  piecewise-quadratic  basis  for  ve¬ 
locity  with  piecewise-linears  for  depth  would  be  a  logical  choice  but  contain 
an  oscillatory  mode  for  velocity  (not  for  depth).  The  nonlinear  terms  in  many 
applications  of  the  shallow  water  equations  for  rivers  and  estuaries  are  quite 
small.  In  our  problem,  however,  the  velocities  are  large  and  the  gradients  steep. 
Thus  the  nonlinearities  are  very  important  and  selecting  the  baisis  to  exclude 
spurious  modes  for  the  linear  equations  does  not  ensure  oscillations  will  not  be 
generated  by  the  nonlinear  convective  terms. 

The  discontinuities  in  h  and  p  at  the  hydraulic  jump  could  be  rep¬ 
resented  exactly  by  means  of  a  discontinuous  finite  element  basis.  In  a  non- 
dissipative  formulation  with  standard  elements  “Gibbs  type”  oscillations  will 
be  generated  at  the  jump.  However,  dissipation  is  present  in  the  physical  prob¬ 
lem  and,  as  noted  previously,  this  effect  is  significant  in  the  neighborhood  of 
the  hydraulic  jump.  Accordingly,  it  is  common  practice  to  include  an  artifi¬ 
cial  dissipation  in  the  governing  equations  or  a  numerical  dissipation  arising 
from  the  discretization  method.  This  added  dissipation  is  active  at  the  hy¬ 
draulic  jump  and  effectively  “smears”  the  jump  over  a  few  mesh  intervals.  The 
higher  order  dissipative  term  then  implies  that  the  upstream  and  downstream 
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boundary  conditions  can  now  be  consistently  treated  and  the  solution  h,  p  to 
the  dissipative  equations  is  continuous.  In  turn,  this  implies  that  we  can  use 
standard  C®  Lagrange  basis  functions  in  the  finite  element  analysis. 


3.4.2  Artificial  Dissipation 


The  form  for  artificial  dissipation  was  constructed  based  on  one¬ 
dimensional  convection  diffusion  models  and  such  that,  throughout  the  su¬ 
percritical  reach,  the  coefficient  closely  matches  the  analytic  wavespeed.  Ac¬ 
cordingly,  we  add  to  (3.1)  -  (3.2)  the  dissipative  operator 


ds 


(3.13) 


where  a  =  ^ji(h)As  in  (3.1)  and  a  —  in  (3.2),  and  is  a  weighting 
parameter.  Adding  the  dissipative  term  (3.13)  to  the  pertubation  equations 
(3.1),  (3.2)  and  integrating  by  parts  in  the  weighted  residual  condition  the 
additional  integrals  in  the  weak  statement  (3.10)  -  (3.11)  are,  respectively. 


and 


juj  dh  dw\ 
[j(h)]  ds  ds  ) 


ds  —  0As 


|u|  dh 
j{h)ds^ 


L7(h)]2  ds  ds 


dp  dw 

—  ~ds-  0As 


1^1  dp 
(j(h)P  ds 


(3.14) 


This  choice  of  dissipation  is  motivated  by  an  analysis  of  the  corre¬ 
sponding  hyperbolic  system  for  the  one-dimensional  spillway  problem.  The 
characteristic  differential  equations  for  the  hyperbolic  system  (3.1)  -  (3.2)  de¬ 
scribing  the  one- dimensional  spillway  problem  are  given  by  (3.9).  The  char¬ 
acteristic  coordinates  are  obtained  as  the  two  families  of  solutions  ^(sj)  = 
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constant  and  Ti{s,t)  =  constant  of  these  characteristic  equations,  in  the  spe¬ 
cial  case  of  «  =  0  the  linearized  shallow  water  equations  are  obtained  and  the 
characteristic  differential  equations  (3.9)  reduce  to 


The  governing  wave  equation  simplifies  in  the  characteristic  frame  and,  in 
the  linear  case  (3.15),  disturbances  propagate  unchanged  along  the  respec¬ 
tive  straight-line  characteristics.  For  the  nonlinear  problem  of  interest  here 
the  characteristics  are  curved.  When  characteristics  of  the  same  family  inter¬ 
sect  to  form  a  caustic,  discontinuities  may  be  generated  in  the  solution.  This 
is  precisely  the  situation  in  the  present  case;  for  supercritical  flow  “positive” 
characteristics  intersect  and  the  difference  between  the  upstream  disturbance 
(which  is  being  convected  downstream)  and  the  downstream  boundary  condi¬ 
tion  is  resolved  through  the  hydraulic  jump.  That  is,  the  specified  tailwater 
in  the  stilling  basin  forces  subcritical  conditions  downstream  while  supercriti¬ 
cal  flow  exists  upstream,  and  the  hydraulic  jump  is  formed  as  a  result  of  this 
interaction. 

If  dissipation  is  introduced  in  the  mathematical  model,  then  the  equa¬ 
tions  are  second-order  and  change  type.  The  upstream  and  downstream  condi¬ 
tions  then  can  be  accommodated  by  a  smooth  solution  to  the  governing  (mod¬ 
ified)  equation.  The  result  of  adding  moderate  local  dissipation  is  that  nu¬ 
merical  oscillations  are  suppressed  in  the  subsequent  approximate  solution  but 
the  sharp  jump  is  Jess  well  resolved.  If  the  local  dissipation  is  small  there  will 
be  only  slight  “smearing”.  This  dissipative  behavior  is,  at  least  qualitatively, 
consistent  with  the  real  physical  process.  This  issue  notwithstanding,  there  are 
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several  other  difficulties  that  arise  when  the  discrete  model  is  considered.  Of 
course,  the  discrete  representation  should  reflect  the  qualitative  behavior  of  the 
continuous  mathematical  problem.  Moreover,  the  approximate  solution  should 
converge  in  the  appropriate  norms  as  mesh  size  and  time  step  are  reduced. 
Our  concern  here  is  to  have  realistic  model  results,  even  near  the  hydraulic 
jump  where  the  effect  of  the  “hard”  downstream  boundary  condition  tends  to 
promote  oscillations  and  spurious  “modes”. 

To  demonstrate  the  nature  of  the  problem,  let  us  turn  to  the  fre¬ 
quently  studied  steady  convection-diffusion  problem.  We  intend  to  clarify  the 
questions  related  to  modes  and  oscillations  associated  with  conflicting  end  con¬ 
ditions,  such  as  those  for  the  hydraulic  jump  in  the  spillway  problem.  The 
model  equation  is 


dc  (Pc  ^ 

—  f-r?  =0  0  <  s  <  1 

ds  ds^ 


(3.16) 


with  c(0)  =  0  and  c(l)  =  1.  Here  e  represents  an  artificial  diffusivity  with 
0  <  t  The  solution  to  (3.16)  is  essentially  zero  in  the  upstream  region 

and  rises  abruptly  through  a  boundary  layer  to  satisfy  c  =  1  at  the  downstream 


Introducing  a  standard  Galerkin  finite  element  method  with  linear 
elements  (or  ecmivalently  using  central  differencing),  at  interior  node  i  we  have 


'c,+i  -  c,  _ 


)  (As)"^""""'  2c. -f-c._,)_ 


(3.17) 


For  the  degenerate  limiting  case  £  — ►  0  the  differential  equation  (3.16)  reduces 
to  dcj dx  =  0  so  c  =  constant  is  the  solution.  In  the  present  case  we  have  c  =  0 
if  the  left  condition  is  to  be  satisfied  and  c  =  1  if  the  right  end  condition  is 
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chosen.  This  degenerate  form  of  (3.16)  is  of  lower  order  and  clearly  is  over- 
specified.  The  corresponding  degenerate  form  of  the  non-dissipative  discrete 
problem  (3.17)  with  £  =  0  presents  an  entirely  different  situation.  We  then 
have  simply  Cj+i  —  Ci_i  =  0.  Setting  Cj  =  p*,  the  solution  to  this  second  order 
difference  equation  satisfies  p*'*’*  —  p'"*  =  0,  that  is,  p  =  ±1.  Hence  the  general 
solution  to  the  discrete  problem  is  c,  =  -f  F(— 1)'  where  constants  A.B  are 
determined  from  the  boundary  conditions.  If  J5  =  0  then  the  discrete  model 
corresponds  to  the  first  order  continuous  problem  and  A  is  given  by  the  end 
conditions.  Even  if  B  is  small,  however,  there  will  be  a  superposed  inter-node 
oscillation  of  amplitude  B.  This  internode  oscillation  of  wavelength  2As  then 
represents  a  spurious  mode  resulting  from  the  second-order  approximation  of 
the  first  order  operator.  The  amplitude  B  of  the  oscillation  will  evidently 
depend  upon  the  “inconsistency”  of  the  boundary  conditions.  In  the  present 
case,  if  c(0)  =  c(l)  then  B  is  zero  and  the  oscillation  will  not  be  forced;  for 
c(0)  =  0,c(l)  =  1  we  see  that  there  is  an  oscillation  of  unit  amplitude  forced 
in  the  numerical  solution. 

Addition  of  artificial  dissipation  returns  the  forms  in  (3.16)  and  (3.17) 
with  t  ^  0  and  a  standard  Fourier  analysis  reveals  that  oscillations  can  be 
eliminated  provided 

c  >  |uoAs  ,  uo  >  0  (3.18) 

which  yields  the  familiar  cell  Reynolds  or  Peclet  condition. 

The  particular  choice  of  artificial  dissipation  employed  in  the  present 
work  can  be  motivated  or  interpreted  physically  as  follows.  First  note  that  the 
coefficient  of  the  artificial  viscosity  is  specified  in  (3.13)  as  the  net  transport 
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rate  u/[j(/i)]^.  This  reduces  to  uo  for  the  standard  shallow  water  equations 
which  then  become 


dh  dp  di^h 

dt'^  ds  ~ 


dp  dh  dp 

m  +  <«*■>  ‘  “"'aJ  +  ^““al 


=  QUO 


ds^ 


(3.19) 


Introducing  a  finite  element  approximation  with  a  linear  baais  the  roots  of  the 
“indicial”  equation  for  the  discrete  difference  form  are 


All  of  these  roots  will  remain  positive  if 


(3.20) 


Q 

As 


> 


(3.21) 


and  numerical  oscillations  excited  by  the  hydraulic  jump  will  be  damped.  More¬ 
over,  if  we  choose  a  to  be  proportional  to  As,  this  condition  will  remain  en¬ 
forced  as  the  mesh  size  As  is  reduced.  The  artificial  dissipation  will  be  reduced 
proportionally  with  As  and  the  resolution  of  the  jump  will  improve.  Since  the 
jump  is  a  major  source  of  oscillation,  the  local  damping  will  improve  the  global 
behavior  of  the  approximate  solution. 


While  this  analysis  has  focused  upon  the  linear  case,  the  nonlinear 
terms  can  also  cause  internode  oscillations.  Aliasing  into  lower  wavelengths 
from  these  internode  oscillations  will  also  occur.  Since  the  artificial  diffusion 
terms  will  preferentially  damp  short  wavelength  oscillations  the  nonlinear  in¬ 
duced  oscillations  should  be  moderated  by  this  mechanism. 
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3.4.3  Finite  Element  Approximation 

Introducing  finite  element  expansions  for  h  and  p  in  the  weak  state¬ 
ments  then  yields  the  discretized  equations 
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for  the  continuity  equation,  and 
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3.5  Model  Results 

Model  simulations  were  conducted  using  bed  topography  from  the 
numerical  and  flume  studies  of  Sivakumaran  (51).  This  bed  form  is  given  by 
Z  =  .20e^" where  X  is  the  horizontal  distance  in  meters  from  the  cen¬ 
ter  of  the  crest  and  Z  is  the  elevation  in  meters.  That  is,  this  spillway  bed 
geometry  has  the  form  of  a  normal  distribution  symmetric  about  z  =  0,  with 
a  crest  elevation  of  0.20m.  In  the  present  numericid  experiments,  values  of 
discharge  are  specified  upstream  and  the  tailwater  elevation  downstream.  The 
flow  capacity  of  the  spillway,  as  indicated  by  the  upstream  water  surface  eleva¬ 
tion  is.  therefore,  calculated  by  the  approximate  model  as  a  steady-state,  from 
which  the  predicted  flow  capacity  can  be  compared  with  Sivakumarsm’s  flume 
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results.  The  same  test  is  repeated  with  the  curvature- related  terms  removed. 
This  yields  the  shallow  water  model  (with  channel  steepness  included),  which 
can  be  used  for  comparison  purposes  to  determine  the  effect  of  curvature. 

One  should  note  that  the  flume  data  of  Sivakumaran  does  not  contain 
a  noticeable  jump  since  the  tailwater  elevation  is  too  low.  However,  in  the 
present  test  a  higher  tailwater  is  specified  to  evaluate  the  suitability  of  our 
numerical  scheme  lo  simulate  the  jump  amd  control  the  associated  oscillations. 
The  location  of  the  jump  can  vary  widely  for  slight  changes  in  velocity  and 
depth.  This  sensitivity  is  reduced  in  real  physical  flows  by  the  presence  of 
friction  which  causes  a  steeper  water  surface  slope.  Friction  in  this  model 
is  included  by  Manning’s  expression  with  friction  factor  0.016  corrected  for 
curvature  [21].  The  results  of  test  calculations  are  shown  in  Figures  3.1  -  3.3 
for  flows  with  an  upstream  discharge  of  0.03599m^/s/m  and  a  downstream 
depth  (tailwater)  0.1m,  using  a  mesh  with  46  linear  elements.  Higher  grid 
resolution  is  required  upstream  of  the  spillway  crest  where  there  are  steep 
variations  in  flow  variables.  Along  the  downstream  side  of  the  spillway  the 
depth  changes  gradually  and  fewer  elements  are  required.  The  semi-discrete 
system  is  integrated  numerically  to  a  steady-state  with  a  fixed  timestep  At  = 
0.05  sec  using  an  implicit  scheme. 

The  depth- averaged  velocity  is  given  in  Figure  3.1.  Both  the  velocity 
and  its  spatial  variation  are  relatively  small  upstream  of  the  crest  which  con¬ 
firms  that  the  effect  of  the  nonlinear  terms  is  not  significant  here.  Downstream 
of  the  crest,  the  velocity  is  large  and  a  small  perturbation  in  depth  imposed 
upon  this  region  will  develop  large  velocity  variations  to  conserve  mass.  If 
dissipation  is  not  sufficiently  large,  the  jump  excites  oscillations  and  a  nonlin- 
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AVERAGE  VELOCmr  PRORLE 
EFFECT  OF  CURVATURE  TERMS 


Figure  3.1:  Average  velocity  profile  —  effect  of  curvature  terms. 

ear  instability  develops.  In  the  calculations  shown  here,  the  coefficient  0  for 
dissipation  was  set  at  0.5  downstream  and  0.25  upstream. 

The  bed  pressure  ratio  profile  is  shown  in  Figure  .3.2.  This  is  the  ratio 
of  the  pressure  associated  with  centrifugal  effects  to  hydrostatic  pressu™  alone. 
Throughout  the  region  in  which  the  bed  is  convex,  —0.24  <  x  <  0.24,  an  uplift 
is  evident,  reaching  a  maximum  of  0.36  a  distance  0.10m  downstream  of  the 
crest.  In  the  concave  portion  of  the  downstream  bed  the  pressure  gradually 
increases  and  reaches  0.55  near  the  “toe"  of  the  spillway.  It  is  apparent  that 
the  pressure  distribution  is  definitely  not  hydrostatic. 

Included  in  Figures  3.1  and  3.3  are  dashed  lines  indicating  the  model 
results  when  the  curvature  is  set  to  be  identically  zero.  This  corresponds  to  the 
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PRESSURE  RATIO  PROFILE 

CENTRIGUGAL/HYOROSTATIC  PRESSURE 


Figure  3.2:  Pressure  ratio  profile  —  centrifugai/hydrostatic  pressure. 

standard  shallow  water  result  for  a  steep  slope.  Of  particular  importance  is  the 
observation  that  the  shallow  water  equations  predict  a  higher  upstream  water 
surface  elevation  to  allow  this  discharge  to  pass  over  the  spillway.  That  is,  the 
shallow  water  equations  predict  a  lower  spillway  capacity  than  the  equations 
including  curvature. 

The  capacity  of  the  spillway  is  directly  related  to  the  energy  required 
to  pass  a  particular  flow  rate.  The  specific  energy,  defined  by  £  s  i  (>^)^  + 
is  the  mechanical  energy  from  the  Bernoulli  expression  referenced  to  the 
datum  of  the  crest  elevation.  Here  (j^)  is  the  surface  velocity  and  is  the 
water  surface  elevation  above  the  crest.  The  specific  energy  indicated  by  the 
flume  was  0.072m.  The  numerical  model  results  were  0.083m  for  the  steep-slope 
standard  shallow  water  equation  and  0.072m  for  the  model  including  curvature. 
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WATER  SURFACE  PROFILE 
EFFECT  OF  CURVATURE  TERMS 
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Figure  3.3:  Water  surface  profile  -  effect  of  curvature  terms. 

Downstream  of  the  crest,  the  calculated  water  surface  for  the  shallow 
water  equations  is  slightly  steeper  than  the  results  with  curvature  included 
and  the  flow  becomes  progressively  more  shallow.  This  is  more  easily  seen  in 
the  velocity  plot,  where  excluding  curvature  produces  velocities  that  are  too 
large.  In  the  concave  region  this  effect  is  reversed  and  the  velocity  predicted 
by  the  steep-slope  shallow  water  equation  drops  below  that  of  the  model  with 
curvature  effects  included.  More  energy  is  needed  to  pass  this  discharge  and  the 
effective  loss  of  energy  downstream  of  the  crest  causes  the  predicted  location  of 
the  jump  to  move  upstream  as  compared  to  tne  results  with  curvature  included. 

The  discrepancy  in  upstream  results  from  the  two  models  is  not  at¬ 
tributable  to  centrifugal  effects  as  they  are  not  significant.  This  behavior  is  due 
to  the  bed-normal  measurement  of  depth.  In  a  concave  region  the  bed-normal 
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directions  at  the  ends  of  differential  segment  converge  as  one  moves  away  from 
the  bed,  so  that  the  differential  volume  per  unit  width  is  smaller  than  would 
be  calculated  by  simply  multiplying  the  depth  and  length  along  the  bed.  The 
volume  (volume  per  unit  width)  calculations  will  be  too  large  in  the  concave  re¬ 
gions  and  small  in  the  convex  regions.  Thus  in  the  concave  regions  the  standard 
equations  require  a  larger  pressure  variation  for  the  flow  acceleration.  Down¬ 
stream  of  the  crest  this  effect  is  smaller  as  the  depth  is  shallow  and  the  water 
surface  is  nearly  parallel  to  the  bed.  Here  the  centrifugal  forces  become  mere 
important.  Results  including  curvature  effects  show  a  slightly  less  steep  water 
surface  slope  as  a  result  of  the  adverse  centrifugal  pressure  gradient.  Down¬ 
stream  of  location  0.4m  the  centrifugal  pressure  gradient  is  positive,  and  the 
results  with  curvature  included  predict  a  flatter  water  surface  slope.  This  in 
turn  moves  the  jump  downstream.  Note  that  both  of  these  model  equations 
are  written  for  bed-fitted  coordinates.  The  St.  Venant  equations  on  the  other 
hand  are  written  for  momentum  and  mass  conservation  horizontally  and  as 
such  offer  a  much  poorer  comparison  on  the  spillway  section  [21]. 

The  capability  of  the  model  to  treat  the  hydraulic  jump  is  also  re¬ 
vealed  in  these  results.  The  effect  of  mesh  refinement  is  examined  in  Figure 
3.4.  A  graded  mesh  with  elements  of  size  0.05m  near  the  jump  was  halved  and 
the  computation  repeated.  Since  the  local  artificial  dissipation  is  proportional 
to  mesh  spacing  the  dissipation  is  reduced  accordingly.  In  both  instances,  the 
jump  occurs  over  three  elements  and,  of  course,  the  steepness  of  the  approx¬ 
imation  to  the  jump  improves  as  the  mesh  is  refined.  There  is  little  spatial 
oscillation. 
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WATER  SURFACE  PRORLE 
EFFECT  OF  REFINEMENT 
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Figure  3.4;  Water  surface  profile  —  effect  of  refinement. 

3.6  Conclusions 

The  finite  element  model  described  here  results  from  a  shallow  water 
perturbation  analysis  including  bed  curvature.  In  the  event  that  the  curvature 
is  assumed  to  be  identically  zero,  the  model  degenerates  to  the  standard  shallow 
water  equations  (without  the  mild-slope  assumption}.  By  computing  the  flow 
over  a  “bump"  composed  of  a  bed  of  continuous  curvature  for  which  flume 
results  were  available  we  observe  that; 

1 .  The  full  equation  model  yields  an  accurate  prediction  of  the  reservoir  ele¬ 
vation  (and  thus  spillway  capacity)  while  the  results  excluding  curvature 
effects  are  much  poorer. 
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2.  The  results  excluding  curvature  effects  compare  quite  closely  with  the 
more  complete  model  along  the  supercritical  spillway  face.  However, 
subtle  differences  are  sufficient  to  significantly  shift  the  predicted  location 
of  the  hydraulic  jump. 

3.  The  resolution  of  the  jump  in  this  model  is  spread  over  2  or  3  elements. 
Thus  the  steepness  of  the  jump  can  be  improved  by  grid  refinement 
(within  the  limits  imposed  by  nonlinear  effects). 

In  the  next  chapter  we  utilize  these  results  to  extend  the  formulation  to  a 
two-dimensional  model. 


Chapter  4 


Model  Extensions 

4.1  Introduction 

The  complexity  of  the  flow  is  more  evident  in  two  dimensions.  The 
lateral  distribution  of  flow  and  oblique  standing  waves  can  be  reproduced.  In 
the  subcritical  case  a  boundary  disturbance  propagates  both  downstream  and 
upstream  and  so  the  shallow  water  equations  cannot  produce  a  standing  wave 
here.  In  the  supercritical  region  a  boundary  disturbance  is  swept  downstream 
forming  a  standing  wave  or  wake.  The  fluid  velocity  is  greater  than  the  speed 
of  the  free-surface  wave  celerity.  Poor  approach  conditions  result  in  standing 
waves  within  the  spillway  apron  section.  The  concentration  of  flow  in  nar¬ 
row  regions  can  cause  circulation  and  subsequent  damage  downstream.  These 
problems  can  be  relieved  with  improved  flow  training  using  well  designed  abut¬ 
ments.  The  flow  features  must  be  modeled  to  assess  the  design,  requiring  a 
model  with  both  longitudinal  and  lateral  resolution. 

This  chapter  is  concerned  with  the  development  of  a  finite  element 
model  that  can  reasonably  address  these  conditions  and  the  nonhydrostatic  in¬ 
fluences.  Section  4.2  details  the  development  of  an  improved  convection  scheme, 
4.3  the  depth-integrated  equations,  4.4  is  the  application  of  the  convection  equa¬ 
tion  to  these  equations,  and  Sections  4.5  and  4.6  the  viscous  stresses  and  bed 
drag,  respectively. 
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4.2  Finite  Element  Treatment  -  Convection  Scheme 


The  focus  of  this  investigation  is  on  the  curvilinear  formulation,  but 
an  adequate  system  to  address  oscillation  control  is  necessary  to  make  this 
calculation.  As  discussed  in  the  previous  chapter  on  the  finite  element  approx¬ 
imation  of  the  one-dimensional  problem,  we  must  be  concerned  with  possible 
oscillations  resulting  from  both  spurious  modes  and  nonlinear  effects.  An  arti¬ 
ficial  diffusion  method  was  successfully  introduced  in  Section  3.4.2,  where  the 
added  dissipation  was  developed  based  on  physical  arguments.  This  system, 
while  sufficient  for  of  the  one-dimensional  equations  in  preliminary  tests,  is  too 
simplistic  for  the  eventual  two-dimensional  model.  We  now  present  a  more 
complex  form  for  two-dimensional  analysis. 

Some  recent  finite  element  approaches  for  the  shallow  water  equations 
include  the  Taylor  weak  form  (Baker  and  lannelli  [4j)  and  the  Petrov  Galerkin 
scheme  (Katopodes  {33])  which,  in  actual  application,  are  almost  identical  (see 
also  [5,  3,  34]),  We  will  illustrate  the  basic  ideas  for  the  Taylor  weedc  form  and 
demonstrate  a  further  improvement  for  our  problem  class. 


Baker's  model  is  based  upon  the  work  of  Donea  [19],  Raymond  and 
Garder  [47]  as  well  as  Lohner,  Morgan  and  Zienkiewicz  [42] .  The  derivation  can 
be  clearly  demonstrated  using  the  one-dimensional  model  transport  equation 


dc  df 


=  0 


(4.1) 


(4.2) 


48 


^m+l 


=  c”'  +  At 


dc”' 

dt 


+ 


At*  d'^c”' 

2  di^ 


(4.3) 


where,  the  superscript  indicates  the  time  step,  i.e.,  t  =  mAt.  Motivated  by 
the  well  known  Lax-WendrofF  approach  in  finite  differences,  we  may  utilize 
relationships  derived  from  (4.1) 


dc  _  df 
dt  ds 


(4.4) 


dPc 

di^ 


—  (  ^ 


(4.5) 


with  0(1,02  parameters.  Following  the  Raymond  and  Garder  scheme  Baker 
and  lanelli  (4]  choose  a  combination  of  (4.5)  and  (4.6)  with  qi  =  02  =  a-  The 
associated  modified  equation  is 


^  ^  (a—  4.  42—^ 

dt  ds  2  ^s  \  dt  ds } 


=  0 


(4.6) 


Katopodes  [33]  arrives  at  essentially  the  same  result  via  the  Petrov 
Galerkin  approach.  The  test  function  is  modified  to  include  derivatives  of  the 
Galerkin  test  function  <t>  as 


~  (i>  +  QA 


d4> 

ds 


(4.7) 


where,  is  the  new  test  function  and  q  is  a  parameter  with  unit  of  1/time. 
Using  (  4.7)  in  the  weighted  residual  statement  for  {  4.1  ), 


/.( 


dc  d<i>  .  d<i>  f dc  dc\\  ,  r  „ 
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where,  is  the  domain  and  dCl  is  the  boundary.  This  is  similar  to  the  form 
proposed  by  Dendy  [18]  and  is  based  upon  the  Streamline  Upwind  Petrov 
Gaderkin  (SUPG)  scheme  of  Hughes  and  Brooks  [29]. 

In  actual  application  of  this  method  to  the  shallow  water  equations 
some  deficiencies  are  apparent.  The  primary  problem,  for  the  hydrodynamic 
conditions  we  wish  to  simulate,  is  that  the  dissipative  mechanism  incorporated 
in  this  method  is  nonexistent  near  critical  flow  conditions;  i.e.,  where  the  flow 
undergoes  transition  from  subcritical  to  supercritical  flow,  or  vice  versa.  This 
allows  oscillations  to  gather  at  the  spillway  crest  or  near  a  jump.  The  standard 
one-dimensional  shallow  water  equations  illustrate  this  point. 


ten 


where. 


These  shallow  water  equations  in  nonconservative  form  may  be  writ¬ 


p  =  uh 


c  =  {gh)^ 


Now,  applying  the  Petrov  Galerkin  approach  to  (4.9)  and  discretizing  we  obtain 
the  corresponding  finite  element  problem.  It  is  instructive  to  examine  the 
linearized  steady  problem  for  A  =  constant,  A.  Then  the  discrete  system  on 
a  patch  at  node  j  becomes 

^A^Q,  -  aAA^^Q,  =  0 


(4.10) 
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where, 

where  we  assume  the  numerical  solution 


h,  =  hp>  (4.11) 

Pj  =  VP^  (4.12) 


j  is  a  nodal  index  and  ( * )  above  a  variable  indicates  the  amplitude.  (These 
operators  are  similar  to  centered  difference  operators  for  first  and  second  deriva¬ 
tives.)  We  may  simplify  the  system  by  defining  P  such  that 

PAP-‘  =  A  (4.13) 

where, 

^  =  f  n  f 
0  A2 

Ai  =  u  -f  c 
A2  =  u  —  c 


Here  Aj  and  A2  are  the  eigenvalues  of  A,  u  and  c  are  the  constant  values  of 
velocity  and  celerity  that  form  A,  and  P  and  P“*  are  made  up  of  the  right  and 
left  eigenvectors  of  A.  We  find 


P  = 


c 


A2 

A| 


(4.14) 


-1  1 

— Aj  A2 


(4.15) 


Equation  (4.10)  may  be  written  as 


-P-^APSq-aP-^APP-'APS^q  =  0 

£i 


or  after  simplification 


-P^Q  -  aAPS^Q  =  0 


(4.16) 


The  equation  set  for  each  patch  will  then  be 


{\\x: 

For  nontrivial  solutions  of  h  and  p  to  exist,  the  determinant  of  the  coefficient 
matrix  must  equal  zero 


X^US-aXiSA  -l6  +  aXiS^ 
Xihs-aX-iSA  -XS-^-aXiS^ 


(4.18) 


from  which  we  have  the  roots 


aAj  +  i  aA2  4-  X 

p  -  1,1, -7 - f,-r - f- 

aAj  -  -  qAj  -  ^ 


(4.19) 


The  value  1  is  a  multiple  root  and  is  the  correct  solution;  i.e.,  the  correct 
solution  for  a  single  boundary  condition  is  that  the  velocity  and  depth  are 
constants;  the  other  values  of  p  are  spurious  numerical  roots.  If  these  spurious 
roots  are  negative,  a  node-to-node  oscillation  may  develop.  The  value  of  a  that 
eliminates  negative  roots  is 


As  As 


(4.20) 


If  IA2I  is  very  small,  as  one  would  expect  near  a  spillway  crest  (assuming  u  >  0), 
a  would  have  to  be  quite  large.  This  is  indeed  what  happens  when  applying  this 
method.  The  term  a  has  units  of  1 /velocity.  Katopodes  suggests  q  =  ^2^’ 
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where  e  is  a  nondimensional  coefficient  (33].  This  again  does  not  scale  the 
weighting  of  the  Petrov  Galerkin  test  function  properly  for  each  characteristic 
and  oscillations  can  occur  in  this  linearized  example. 


Instead  we  propose  that  the  test  function  be  modified  so  that  it  is 
scaled  by  each  characteristic  magnitude.  We  define  the  new  test  function  as 


where, 


-  dd> 
+  oA— 
as 


(4.21) 


A  = 

A  = 

The  result  for  the  discrete  model  on  a  patch  then  is  approximately 

AP^Q  -  ap-‘APP-‘AP«^Q  =  0 


P‘‘AP 

Ai 
|Ai| 

0 


0 

IAjI  J 


and  to  prevent  oscillations 


(4.22) 


This  method  provides  an  upwinded  Petrov  Galerkin  test  function  that  is  scaled 
for  each  characteristic  and  can  provide  oscillation  control  at  all  Froude  num¬ 
bers,  We  will  substitute  a  =  cAs  , where  c  is  a  nondimensional  parameter, 
in  the  model  so  that  we  need  not  input  values  dependent  on  element  lengths. 
This  method  is  a  variation  on  the  theme  proposed  by  Courant,  Isaacson,  and 
Rees  [15]  for  one-sided  finite  differences.  These  ideas  were  expanded  to  more 
general  problems  by  Moretti  [45],  Chakravarthy  [11],  and  Gabutti  [25]  as  split- 
coefficient  matrix  methods  and  by  the  generalized  flux  vector  splitting  proposed 
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by  Steger  and  Warming  [55].  These  one-sided  differences  based  upon  character¬ 
istic  directions  provide  dissipation  and  therefore  stability.  The  incorporation 
here  is  more  subtle  using  a  modification  of  the  test  function;  also  the  degree  to 
which  the  dissipation  is  apparent  is  controlled  via  the  parameter  c  instead  of  a 
totally  one-side  difference. 


This  approach  is  now  extended  to  the  two-dimensional  case.  We  will 
use  the  condition  of  zero  curvature  in  deriving  our  test  function.  This  is  sig¬ 
nificantly  simpler  than  including  curvature  and  can  be  justified  on  physical 


grounds.  The  shallow  water  equations  (with  and  C2  not  identically  equal  to 
1)  are 


ut  USi  U$2 


(4.23) 


where, 


H  represents  the  additional  terms  which  are  not  associated  with  the  material 
derivative  or  pressure,  and  gj  is  the  gravitational  component  in  the  bed-normal 
direction. 
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The  test  function  we  propose  is 


ij)  - 


(4.24) 


where, 


A 

A  = 
A  = 

Ai  = 

A2  = 
A3  = 

A  = 

p  = 


,  .  d<f>  * 
c  (  — A  +  As? 

</Si 


p-*Ap 

Ai  0  0 

0  Aj  0 
0  0  A3 

u 

z 

u  +  c 

~z~ 

u  —  c 


£•) 


lb 


0 

0 


0 

Al 

|Aal 

0 


p->  =  ^ 


1_ 

2c 


0 
0 

lAsI  J 

-u  0  1 
-Cl -^3  1  0 
-C1A2  1  0 
0  1  -1 
0  C1A2  —Cl  A3 
2c  V  —V 


As  in  the  one-dimensional  case,  Ai,A2,  and  A3  are  the  eigenvalues  of  A,  and  P 
and  P"*  are  made  up  of  the  right  and  left  eigenvectors. 

Similarly,  in  the  S2  direction  we  have: 


B  =  R-*rR 

■  7,  0  0 

r  =  0  72  0 

.0  0  73  J 
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71  = 

72  = 

73  = 

r  = 

R  = 

R-‘  = 


The  terms  F,  R,  and  R“‘  in  the  $2  direction  are  analogous  to  the  relationship 
of  /1,P,  and  P~*  in  the  si  direction. 


4.3  Model  Equations 

The  depth-integrated  version  of  the  momentum  equations  (2.58)  is  of 
some  advantage  in  handling  non-smooth  conditions.  The  weak  form  solution 
for  the  no-curvature  condition  that  one  would  encounter  downstream  of  the 
spillway  can  be  made  to  properly  conserve  momentum  and  mass  through  the 
jump.  Other  forms  of  the  equations  will  not.  In  the  case  in  which  there  is  bed 
curvature,  as  we  discussed  in  Section  3.3,  these  equations  will  contain  additional 
terms  due  to  the  bed  curvature  which  while  finite  through  the  jump  will  make  an 
additional  contribution  that  can  cause  an  error  in  the  jump  location.  Therefore 
in  the  vicinity  of  the  jump  these  equations,  which  properly  conserve  meiss  and 
momentum  for  the  no-curvature  case,  will  only  precisely  conserve  m<iss  in  the 
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curved  bed  state.  Generally,  in  practical  cases  the  strong  jump  is  restricted  to 
the  region  downstream  of  the  spillway  face  where  the  bed  contains  no  curvature. 
These  equations  then  are 


L(Q)  =  Mt  + 


C1C2L  dsi 


5(C2Na)  a(CiNb) 


ds2 


+  7-rH  =  0  (4.25) 

C1C2 


where, 

Q 
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/2(53) 
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and  the  nonhomogeneous  term 
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where. 
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"1  =0 

+  ^  [“S®3  +  ^ 

+53{-|?[/‘ax-a.)  +  ^|ff}  +  CiC25i«7  ^ 

Hs  =  |?T{-i7 

+2^{~l&  [(awF'M  [H/lfe)  ~H} 

+  fT  [“il^  +  ^  [1?^  “  C2|7f^] 

+«2^  {[&  (V)l  ^  +  2Cif  |^<>11 

+53  {-§^7  [^^1  -  ^4]  +  ^1^}  +  CiC252<>7 
With  the  coefficients  a^  ,  6„  and  c,  resulting  from  depth  integration: 
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Kjh  IV  K2 

For  the  coefficients  b^  simply  swap  the  subscripts  1  and  2  in  the  kj  and  /C2  and 
fi  and  /2  terms  of  a; . 


4.4  Application  to  Generalized  Equations 

The  Petrov  Galerkin  formulation  incorporates  a  combination  of  the 
Galerkin  test  function  and  a  non-Galerkin  component  to  control  oscillations 
due  to  convection.  The  test  function  is  that  of  equation  (4.24)  for  the  zero- 
curvature  Ccise,  and  involves  bed  velocities  u  and  v. 


The  weighted  residual  statement  becomes; 


<^Mt  -  -r-;g— N«  -  dQ 

(.l^Si  (,2<)S2  (,lC,2 


+  /  (N«ni  +  Nb^a)  =  0 
JBO 
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(4.26) 

where,  ni  and  02  are  the  components  of  the  outward  normal  vector  to  the 
boundary. 

The  perturbation  partial  differential  system  defines  a  hyperbolic  ini¬ 
tial  boundary  value  problem.  We  determine  the  appropriate  boundary  condi¬ 
tions  relying  upon  the  approach  of  Daubert  and  Graffe  [17]  and  the  discussions 
of  Drolet  and  Gray  {23]  and  Verboom  et.  al.  [58].  Daubert  and  Graffe  use 
the  method  of  characteristics  for  this  determination.  The  theory  shows  that 
the  number  of  boundary  conditions  is  equal  to  the  number  of  characteristic 
half  planes  that  originate  exterior  to  the  domun  and  which  enter  it.  There 
are  two  families  of  characteristic  surfaces  at  a  point  (.Sio,  .820,  ^o)-  The  first  is  a 
cone  which  slopes  in  the  direction  of  flow  and  has  a  radius  of  the  wave  celerity 
multiplied  by  time.  This  generates  the  ring  formed  on  the  free-surface  by  a 
disturbance.  The  second  is  a  plane  which  intersects  the  axis  of  the  cone  and 
represents  the  flow  velocity.  If  a  characteristic  plane  outside  our  domain  inter¬ 
sects  the  boundary  then  the  flow  field  inside  the  domain  is  influenced  by  the 
outside  information  which  we  must  provide  as  a  boundary  condition.  On  the 
other  hand,  if  the  characteristic  plane  intersects  the  boundary  from  inside  the 
domain,  this  is  information  leaving  the  domain  and  no  boundary  condition  is 
needed.  Physically,  the  first  family  is  tracking  a  free-surface  disturbance  and 
the  second  is  tracking  a  fluid  element.  Table  4.1  relates  the  number  of  bound¬ 
ary  conditions  to  the  direction  of  flow  for  sub-  or  supercritical  flow  conditions. 
Ur,  is  the  velocity  component  normal  to  the  boundary. 

The  boundary  conditions  we  implement  in  this  model  are: 
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Table  4.1:  Number  of  boundary  conditions  required. 


Flow  direction 

Subcritical  Flow 

Supercritical  Flow 

tin  <  0  (inward  flow) 

2 

3 

u„  >  0  (outward  flow) 

1 

0 

u„  =  0  (no  flow) 

1 

1 

•  A  slip  flow  boundary  at  a  solid  wall  is  common  in  shallow  water  ap¬ 
plications.  This  corresponds  to  the  condition  u„  =  0  in  Table  4.1.  The 
boundary  condition  in  the  model  is  imposed  through  the  weak  statement. 
The  implementation  is  as  follows; 

Mass  Conservation  Equation 

/  <^(pcin,  +  gcanj)^/  =  0  (4.27) 

Jan 

si  Momentum  Equation 

+  =  0  (4.28) 

S2  Momentum  Equation 

•  The  upstream  boundary  condition  is  an  essential  boundary.  If  the  flow 
here  is  subcritical  then  the  flow  components  are  specified  as  either  u  and 
V  OT  p  and  q.  If  the  flow  is  supercritical  the  depth  is  also  specified. 

•  The  downstream  boundary  condition  is  unspecified  if  the  flow  is  supercrit¬ 
ical  there.  If  the  flow  is  subcritical  the  depth  is  specified  by  substituting 
the  specified  depth  in  the  boundary  integral  terms  of  the  momentum 
equations. 
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4.5  Viscous  Stress 


The  viscous  stress  is  of  less  importance  in  the  calculation  of  the  free- 
surface  except  as  the  flow  passes  through  the  hydraulic  jump.  The  contribution 
is  small  when  the  flow  is  smooth.  We  base  our  stress  calculation  upon  the 
velocity  and  metric  gradients  at  the  bed,  in  this  way  we  do  not  produce  a 
lateral  or  longitudinal  stress  by  our  assumed  vertical  velocity  profile.  The 
stresses  (see  Boresi  and  Chong  [8]  }  then  are: 

-  = 

(  \  du  \  dv  V  8(2  “  ^Ci  \  /A  01  \ 

<^2i,<7i2  —  I  7"  Q 7“ a - 7T“a - — )  (4-31) 

\C20s2  osi  CiCa^^i  C1C2052/ 


0-22  =  2up 


dv  u 

^■*2  C1C2 


1C2 


(4.32) 


The  first  subscript  indicates  the  face  upon  which  the  stress  acts  and  the  second 
is  its  direction.  When  applied  to  our  numerical  system  we  obtain  the  following 
additional  contributions; 


si  Momentum  Equation  Contribution: 


Jn  p  VCi  dsj 


1 


(4.33) 


S2  Momentum  Equation  Contribution: 


These  contributions  are  made  to  the  left  hand  side  of  equation  (4.26).  The 
stress  jump  across  interelement  boundaries  is  set  to  zero  and  along  the  domain 
edge  a  stress-free  boundary  is  assumed. 
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As  an  estimate  of  these  turbulent  eddy  viscosities  generated  by  the 
bottom  friction,  we  use  an  empirical  formula  (49,  12) 

=  Cb/i/’  i«i  (4.35) 

where,  is  the  turbulent  viscosity,  Cb  is  a  coefficient  that  varies  between  0.1 
and  1.0,  and  /  is  the  Darcy- Weisbach  friction  coefficient  (see  e.g.  King  and 
Brater  [36]).  The  turbulent  contribution  to  viscosity  is  generally  much  larger 
than  the  molecular  viscosity. 


4.6  Bed  Drag 

We  adopt  the  common  practise  in  hydraulic  engineering  of  using  an 
empirical  relationship  developed  by  Manning  and  extended  to  curved  beds  in 
one  dimension  by  Dressier  and  Yevjevich  [21].  We  shall  extend  this  to  two 
dimensions.  The  original  approach  uses  the  Chezy  coefficient,  Ch,  defined  for 
steady  flow  and  for  determination  of  the  average  flow  over  a  cross  section.  In 
the  “fully  rough”  regime,  the  bed  stress  may  be  defined  by  tjbfi  —  ~\pu^,  where 
A  and  C/,  are  functions  of  the  size  of  the  flow  channel  and  the  bed  roughness 
and  independent  of  Reynolds  number  and  so  the  Manning’s  relationship  is 
applicable.  The  term  crbtJ  will  be  the  bed  drag  force  addition  to  equation 
(4.26).  Here  we  apply  Manning’s  relationship  of  A  =  where  R  is  the 
hydraulic  radius  (ratio  of  cross-sectional  area  to  wetted  perimeter)  and  n  is 
Manning’s  coefficient. 

We  shall  demonstrate  the  derivation  of  the  force  terms  for  the  sj 
momentum  equation.  Consider  the  infinitesimal  volume  in  Figure  4.1.  The 
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Figure  4.1:  Infinitesimal  volume  for  calculation  of  bed  drag  force. 


bed  surface  area  is 

A  =  <’i<’2A5iAJ2 

so  the  total  bed  drag  force  is  given  by 

fsl  =  +  V*)*  CiC2'^S,AS2 


The  corresponding  volume  is 


y  =  f  jl  (S3)  J2  (S3)  </s3CiC2ASiAs2  =  a7ClC2ASi  As2 
Jo 

and  the  force  per  unit  volume  then  is 

/»!  9^  -  { ‘2  ,  *2^2 

T7  =  -'"sr"  (“  +  *’ ) 


V  -  '■Jli 

The  hydraulic  radius  may  be  approximated  by 


(4.36) 


(4.37) 


(4.38) 


(4.39) 


^  V 
ft  SS  —  =  07 
A 


(4.40) 


Therefore,  the  force  per  unit  mass  as  it  appears  on  the  right  hand  side  of  the 
S)  momentum  equation  becomes: 
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Fbt  =  - 


gnV((p^-4-(?^)^ 

hoj 


and  similarly  for  the  S2  momentum  equation: 


Fb2  =  - 


gn^q{ip^  +  q^)'^ 
hi4 


(4.41) 


(4.42) 


Once  again,  these  terms  will  appear  as  additions  in  (4.26)  for  the  bed  drag 
force. 


4.7  Model  Description 


The  finite  element  approximation  for  equation  (4.26  )  becomes: 


i  (N.m  +  Npn,)  dl] 


^  CiG 


(4.43) 

0 


where,  the  symbol  *  indicates  the  discrete  value  of  the  quantity  and  the  sub¬ 
script  indicates  a  particular  test  function.  The  geometry  and  flow  variables  are 
represented  using  the  finite  element  basis;  e.g.. 


Q  =  E 

3 

and  the  function  is  defined  as 


Vi 


e 


ddii  ; 

—A  +  As2 

(isi 


We  use  quadrilateral  bilinear  elements  with  nodes  at  the  element  cor¬ 
ners,  the  local  element  coordinates  are  shown  (see  e.g.  [6])  in  Figure  4.2.  The 
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Figure  4.2:  Local  bilinear  element. 


grid  intervals  Asi  and  As2  in  the  physical  domain  are  chosen  in  the  same 
manner  as  Katopodes  [33]: 


Aso  =  2 


9  j.  (2h.\ 

[Ud  \dv) 

■  -  (l)'*(S) 


(4.44) 


(4.45) 


From  Taylor  series  arguments,  the  temporal  derivative  may  be  repre¬ 


sented  as 


^  At  At  dt  ^  ^ 


(4.46) 


where,  the  subscript  j  indicates  a  particular  node  location,  the  superscript 
indicates  the  time  step,  and  Tj  indicates  the  truncation  error. 


- 1)  -{■’-Dw 
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If  7  =  2  defines  a  second-order  backward  difference  in  time;  if  7  =  1  the 
standard  first-order  backward  difference  is  obtained  (both  schemes  are  uncon¬ 
ditionally  stable  for  the  linear  problem). 

We  address  the  nonlinearities  in  equation  (4.44)  by  using  the  Newton 
Raphson  technique  (see  [9]).  The  residual  for  a  particular  test  function  4>i 
in  each  set  of  equations  is  forced  toward  zero  by  the  iteration:  for  /:  =  1, 2, . . . 
solve  the  Jacobian  system 

=  -R*  (4.48) 

=  Q;  +  AQj  (4.49) 

where  k  indicates  the  iteration,  i  is  the  test  function  index  and  j  is  the  nodal 
index. 

Equation  (4.48)  represents  a  system  of  linear  algebraic  equations  that 
must  be  solved  for  each  iteration  and  time  step.  A  “Profile”  solver  incorporates 
efficient  coefficient  matrix  storage  and  is  implemented  in  the  present  study.  In 
this  method  the  upper  triangular  portion  of  the  coefficient  matrix  is  stored  by 
columns  and  the  lower  by  rows.  The  zeros  outside  the  profile  are  not  stored  or 
involved  in  computation.  The  necessary  arrays  are  then  a  vector  comprised  of 
the  columns  of  the  upper  triangular  portion  of  the  coefficient  matrix,  another 
for  the  rows  of  the  lower  portion,  and  a  pointer  vector  to  locate  the  diago¬ 
nal  entries.  Triangular  decomposition  of  the  coefficient  matrix  is  used  in  a 
direct  solution.  The  program  to  construct  the  triangular  decomposition  of  the 
coefficient  matrix  uses  a  compact  Crout  variation  of  Gauss  elimination. 


A  r.k 
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Chapter  5 


Validation  and  Practical  Simulation 


Numerical  tests  are  conducted  in  two  stages.  In  Section  5.1  we  make 
validation  tests  against  several  standard  benchmark  cases.  Then  in  Section 
5.2  we  demonstrate  the  model’s  capability  to  predict  actual  flow  conditions. 
Here  we  make  comparisons  to  flume  results  d  compare  these  results  to  those 
predicted  using  other  shallow  water  equation  models  of  approximately  the  same 
computational  complexity. 

5.1  Model  Validation 

Several  simple  cases  with  analytic  solutions  are  first  considered  to 
validate  the  scheme  and  program.  The  basis  of  these  is  the  Bernoulli  Equation 
written  along  the  water  surface.  For  smooth  solutions  with  no  dissipation, 
energy  along  a  streamline  should  be  constant.  The  Bernoulli  equation  along 
the  water  surface  may  be  written: 

where. 


Em  =  mechanical  energy  (in  units  of  length,  e.g.  ft-lbs  per  lb  of  water) 

z,  =  water  surface  elevation 

V,  =  velocity  magnitude  at  the  water  surface 
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In  testing  the  precision  of  the  model  we  define  the  error  in  terms  of  deviation 
from  a  constant  mechanical  energy  as  follows; 


Bn. 

l|en.ll 


loBmjl  {h)j2(h)d3ids2 

/nil  ih)j2ih)dsxds2 

fo  (Bm  -  Et^y  jiih)j2(h)dsids2 
/nii(*)ii(^)‘^«i<^^2 


(5.2) 

(5.3) 


So  that  Em  is  the  average  mechanical  energy  and  ||£mll  represents  an  error  in 
terms  of  mechanical  energy. 


5.1.1  Constant  Curvature. 

The  first  test  is  for  a  constant  curvature.  These  tests  are  conducted 
for  constant  curvature  in  the  Sj  direction  with  no  curvature  in  the  sj  direction, 
and  then  for  constant  curvature  in  the  $2  direction  with  no  curvature  in  the  sj 
direction.  All  tests  are  with  straight  walls  for  both  supercritical  and  subcritical 
conditions. 

The  radius  of  curvature  chosen  is  30m,  and  the  channel  length  is  20m 
divided  into  elements  of  length  Im.  The  width  is  4m  divided  into  Im  wide 
elements.  In  addition  to  the  error  in  mechanical  energy,  we  plot  the  depth  and 
bed  velocity  along  the  channel  centerline  compared  to  the  analytic  solution  for 
the  calculated  average  energy.  We  show  these  plots  for  flow  in  the  Si  direction; 
the  $2  plots  were  indistinguishable  from  the  S)  results. 

The  subcritical  test  conditions  are  shown  in  Table  5.1.  For  subcritical 
flow  two  upstream  boundary  conditions  and  one  downstreaim  boundary  condi¬ 
tion  are  specified.  Typically  p  and  q  are  given  upstream  and  water  surface 
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Table  5.1:  Circular  channel,  subcritical  flow. 


Input  Parameters 

9 

mjhhwnfimm 

n 

0 

I't 

0 

At 

2sec 

Simulation  Time 

lOOsec 

Iterations 

2 

c 

0.01 

Boundary  Conditions 

Upstream 

P 

1.048  m^/s/m 

9 

0 

Downstream 

Tailwater  Elevation 

1.0439m 

Table  5.2;  Circular  channel  subcritical  flow  results. 


Ikmil 

■Sl 

1.046120 

4.024x10-“ 

S2 

1.046122 

4.084x10-*^ 

elevation  downstream.  The  results  of  the  model  run  at  the  end  of  the  simula¬ 
tion  are  shown  in  Table  5.2.  The  mechanical  energy  in  the  model  is  very  close 
to  a  constant.  For  subcritical  flow  the  energy  should  be  close  to  the  variation 
in  elevation,  so  it  is  apparent  that  the  error  in  elevation  is  less  than  0.001m. 
A  plot  of  the  centerline  profiles  of  bed  velocity  and  water  surface  is  shown 
in  Figure  5.1.  The  water  surface  dips  over  the  crest  as  one  would  expect  for 
subcritical  flow  and,  of  course,  the  velocity  is  a  maocimum  there.  A  comparison 


ELEVATION  .  M  AND  VELOCITY  .  M/S 


70 


CENTERUNEPRORLE 


CIRCULAR  BED 


Figure  5.1:  Circular  channel,  subcritical  flow;  bed  velocity  and  water  surface 
elevation. 
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CENTERUNE  DEPTH  COMPARISON 
CIRCULAR  BED 


Figure  5.2;  Circular  channel,  subcritical  flow;  comparison  of  depth  from  model 
with  analytic  solution. 

between  the  numerical  model  and  analytic  results  is  shown  in  the  two  subse¬ 
quent  figures.  Figure  5.2  provides  a  comparison  of  depths  and  Figure  5.3  is  for 
bed  velocity.  The  comparison  is  extremely  close,  and  since  the  energy  error  is 
small  we  know  the  comparison  over  width  is  also  good. 

The  test  for  supercritical  flow  considers  the  treatment  and  effect  of 
the  downstream  boundary.  Since  the  flow  is  supercritical  at  this  boundary  no 
condition  is  specified.  The  upper  portion  is  subcritical  so  that  only  p  and  q 
must  be  specified  upstream.  The  flow  changes  to  supercritical  near  the  crest. 
The  water  surface  rises  or  falls  from  initial  data  to  the  level  at  which  the 


steady-state  discharge  can  be  maintained. 

The  downstream  depth  is  c^Jculated  by  the  model  and  as  a  result  is 


CENTERUNE  VELOCITY  COMPARISON 
CIRCULAR  BED 


10.  15.  20. 
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Figure  5.3:  Circular  channel,  subcritical  flow;  comparison  of  the  bed  velocity 
from  the  model  with  analytic  solution. 
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Table  5.3:  Circular  channel,  supercritical  flow. 


Input  Parameters 

9 

9.80m/sec^ 

n 

0 

>'t 

0 

At 

Isec 

Simulation  Time 

lOOsec 

Iterations 

O 

im 

t 

0.01 

Boundary  Conditions 

Upstream 

P 

1.048m^/s/m 

9 

0 

Table  5.4:  Circular  channel  supercritical  flow  results. 


Er. 

Ikmil 

0.6973918 

1.911x10-^ 

S2 

0.6973964 

1.912x10-^ 

somewhat  less  stable  so  the  time  step  is  reduced  accordingly.  Table  5.3  shows 
the  input  parameters  for  the  calculation. 

The  error  compared  with  Bernoulli’s  equation  is  shown  in  Table  5.4. 
As  expected  the  error  is  greater  than  for  the  subcritical  case  but  it  is  still 
relatively  small. 

The  actual  flow  centerline  profiles  are  shown  in  Figure  5.4.  The  water 
depth  decreases  until  it  becomes  supercritical  near  the  crest.  The  bed  velocity 
is  nearly  7m/s  at  the  downstream  boundary.  The  comparison  to  the  analytic 
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Figure  5.4:  Circular  channel,  supercritical  flow;  bed  velocity  and  water  surface 
elevation  profiles. 


result  is  shown  in  Figures  5.5-5.6.  Once  again,  the  model  yields  results  very 
close  to  the  analytic  solution. 


5.1.2  Variable  Curvature. 

For  a  more  general  test  we  consider  a  variable  curvature  case.  Specif¬ 
ically,  we  return  to  the  flume  geometry  of  Sivakumaran  [51,  52,  53]  given  by 
Z  =  0.20cI-3<  ^^*1,  where  x  is  the  horizontal  distance  and  Z  is  the  elevation. 
The  test  input  parameters  are  shown  in  Table  5.5.  The  simulated  flume  con¬ 
tains  245  nodes  and  192  elements.  The  lateral  resolution  consists  of  4  elements 
of  width  0.05m.  Longitudinally,  the  resolution  is  concentrated  just  upstream  of 
the  crest,  where  the  greatest  variation  in  depth  occurs.  The  minimum  element 
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Figure  5.5;  Circular  channel,  supercritical  flow:  comparison  of  depth  from 
model  with  analytic  solution. 


Table  5.5:  Variable  curvature  flume. 


Input  Parameters 

9 

9.80m/ sec^ 

n 

C 

I't 

0 

At 

O.OSsec 

Simulation  Time 

20sec 

Iterations 

3 

f 

0.02 

Boundary  Conditions 

Upstream 

P 

0.03599mVs/m 

9 

0 

BED  LONGITUDINAL  VELOCITY  .  M/S 
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CENTERUNE  VELOCITY  COMPARISON 
CIRCULAR  BED 


Figure  5.6;  Circular  channel,  supercritical  flow;  comparison  of  the  bed  velocity 
from  the  model  with  analytic  solution. 
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Table  5.6;  Variable  curvature  flow  results. 


lien.!! 

Si 

0.273170 

2.959x10-^ 

S2 

0.273090 

3.990x10-* 

length  is  0.02m  and  the  maximum  is  0.05m.  With  a  0.05  second  timestep  a 
free-surface  wave  can  traverse  about  3  or  4  of  the  smallest  elements  in  a  single 
timestep;  i.e.,  the  equivalent  CFL  number  [14]  is  about  3  or  4.  Our  initial  solu¬ 
tion  is  a  fairly  poor  description  of  the  steady  solution  so  the  starting  timestep 
is  0.02  seconds.  After  2  seconds  of  simulation  it  was  increased  to  0.05.  The 
results  of  the  simulation  are  shown  in  Table  5.6. 

The  ratio  ||£„||/£„  for  the  variable  curvature  case  is  between  0.00108 
and  0.00146  compared  to  the  supercritical  constant  curvature  result  of  about 
0.00274  and  subcritical  result  of  0.000390.  Thus,  the  variable  curvature  is  no 
worse  than  for  supercritical  flow  over  a  constant  curvature.  The  significant 
difference  is  in  comparison  to  subcritical  flow.  Supercritical  flow  by  nature 
is  less  regular  and  the  results  with  the  model  reflect  this.  These  validation 
studies  confirm  that  the  model  accurately  approximates  the  stated  differential 
equations.  We  now  apply  the  model  to  more  realistic  practical  flow  studies. 

5.2  Comparison  to  Physical  Measurements  in  Flumes 

The  first  test  contains  no  bed  curvature,  but  is  a  supercritical  transi¬ 
tion  in  which  the  disturbances  at  the  wall  propagate  across  the  flume.  This  is 
a  good  test  of  the  numerical  scheme  and  boundary  condition  implementation. 
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Figure  5.7;  Flume  geometry;  supercritical  transition. 

We  next  compare  results  for  our  generalized  set  of  shallow  water  equations,  the 
standard  steep-slope  shallow  water  equations,  and  the  St.  Venant  equations  of 
a  spillway  form  against  data  collected  by  Sivakumaran  [51].  These  methods  are 
all  of  about  the  same  computational  effort.  The  data  set  includes  bed  pressures 
and  water  surface  elevations.  The  final  test  is  of  an  outletworks  flume  that  is 
supercritical  throughout  containing  bed  curvature  and  lateral  transitions.  This 
is  the  most  general  case  that  we  will  test.  In  all  comparisons,  model  resolution 
is  refined  until  no  significant  change  results  from  additional  refinement. 

5.2.1  Supercritical  Transition. 

The  flume  data,  for  this  test,  are  reported  in  Ippen  and  Dawson  [31]. 
The  flume  narrows  from  2ft  to  1ft  wide  using  two  equal  radius  circular  arcs 
eis  shown  in  Figure  5.7.  The  numerical  representation  of  this  flume  extends 


Figure  5.8:  Computational  grid,  supercriticad  transition. 

upstream  20ft  from  the  start  of  the  curve  and  downstream  below  the  transition 
by  about  30ft.  The  extension  upstream  is  to  allow  the  model  to  reach  uniform 
flow  before  the  transition  and  downstream  to  make  sure  the  boundary  condition 
does  not  influence  the  calculations  in  the  area  of  interest.  These  long  extensions 
are  actually  overly  cautious.  The  numerical  grid  is  composed  of  4585  nodes 
and  4314  elements  with  a  maximum  length  of  1ft  and  a  minimum  of  0.06ft. 
Laterally  the  channel  is  broken  into  20  elements  in  the  area  of  interest  and  near 
the  upstream  and  downstream  boundaries  only  6  elements.  The  numerical  grid 
in  the  transition  area  is  shown  in  Figure  5.8.  Other  model  input  parameters 
are  shown  in  Table  5.7.  The  bed  slope  and  roughness  are  chosen  to  provide 
uniform  flow  approaching  the  transition.  The  flow  conditions  are  more  difficult 
than  in  the  previous  test  and  the  value  of  e  was  increased  to  0.10  to  maintain 
stability.  The  value  of  i/,  is  chosen  to  correspond  to  Cb  —  0.4,  see  equation 
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Table  5.7:  Supercritical  transition. 


Input  Parameters 

9 

32.208ft /sec^ 

n 

0.005 

t't 

0.02ft^/sec 

At 

0.03sec 

Simulation  Time 

500sec 

Iterations 

2 

e 

0.10 

bed  slope 

0.0125 

Boundary  Conditions 

Upstream 

P 

0.7161ftVs/ft 

9 

0 

h 

0.0998ft 

(4.35).  The  amplitude  of  the  standing  wave  predicted  by  the  numerical  model 
is  fairly  sensitive  to  Cb  as  one  might  expect. 

A  perspective  view  of  the  computed  water  surface  is  shown  in  Figure 
5.9.  The  transition  causes  a  disturbance  that  reflects  down  the  channel  forming 
a  diamond-shaped  wave  pattern.  A  comparison  of  computed  and  experimental 
water  surface  contours  is  shown  in  Figure  5.10.  The  numerical  model  certainly 
captures  the  overall  features  of  the  flume.  Not  surprisingly  it  is  symmetric 
unlike  the  flume  (since  it  is  difficult  to  control  the  inflow  into  a  physical  flume). 
The  numerical  model  predicts  the  location  of  the  initial  peak  upstream  of  the 
flume’s  peak  by  about  0.5ft  (about  15  %  of  the  transition  length)  and  each 
subsequent  peak  is  increasingly  off.  The  distance  between  the  first  two  peaks  is 
3.5ft  but  the  model  produces  only  2.8ft  separation.  Further  down  the  channel 
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Figure  5.9:  Supercritical  transition:  water  surface,  3-(l  perspective  view. 
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the  wavelength  in  the  flume  is  3.3ft  and  the  numerical  model  gives  3.0ft.  This 
is  a  direct  result  of  the  shallow  water  assumption.  The  various  wavelengths 
actually  should  travel  at  different  speeds  with  the  shorter  ones  traveling  more 
slowly.  In  the  shallow  water  model  all  waves  travel  at  the  speed  of  an  infinitely 
long  wave.  For  example,  from  equation  (1.1)  we  see  that  the  long  wave  celerity 
for  a  depth  of  0.3ft  is  3.10ft/sec.  For  the  longest  wavelength  in  the  flume,  3.3ft, 
the  celerity  is  2.95ft/sec;  the  shorter  wavelengths  are  even  more  significantly 
reduced.  This  higher  wave  celerity  means  that  the  waves  which  originate  at 
a  sidewall  boundary  will  form  a  standing  wave  that  is  swept  less  backwards 
by  the  current.  Overall  though,  within  the  limitations  of  the  shallow  water 
assumptions,  the  comparison  is  reasonable  with  symmetry  preserved  and  the 
shapes  of  the  oblique  standing  waves  demonstrated. 

5.2.2  Spillway  Form 

At  this  point  we  revisit  questions  concerning  the  generalized  shallow 
water  equations  as  differentiated  from  the  more  frequently  applied  St.  Venant 
equations  and  the  generalized  shallow  water  equations  with  curvature  set  to 
zero  (the  standard  steep-slope  shallow  water  equations).  We  do  this  by  recon¬ 
sidering  Sivakumaran’s  flume  results  with  the  improved  model,  continuing  the 
preliminary  study  reported  in  Chapter  3.  Here  we  first  address  longitudinal 
issues  since  the  flume  is  essentially  one-dimensional,  leaving  lateral  and  more 
general  two-dimensional  issues  to  the  next  section. 

While  the  generalized  shallow  water  equations  are  much  more  com¬ 
plicated  than  the  St.  Venant  equations,  we  are  using  bilinear  elements  with  the 
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identical  nunnber  of  gauss  points  for  each  approach.  Hence  the  models  are  of 
approximately  the  same  computational  effort.  If  one  moved  to  a  higher-order 
perturbation  expansion,  much  more  computational  effort  would  be  necessary. 
Previous  studies  by  Dressier  and  Yevjevich  [22]  compare  Dressler’s  model  with 
the  St.  Venant  equations.  Since  the  “Dressier”  equations  gave  markedly  bet¬ 
ter  results,  they  concluded  that  “curvature”  is  quite  important.  However,  the 
St.  Venant  equations  enforce  another  simplification  that  may  contribute  to  the 
difference  -  the  implementation  of  the  mild-slope  assumption.  We  shall  inves¬ 
tigate  the  nature  of  the  discrepancies  between  the  generalized  shallow  water 
equations  and  the  St.  Venant  equations.  The  standard  steep-slope  shallow 
water  equations  will  serve  as  a  means  of  evaluating  this  additional  assumption 
in  the  St.  Venant  form. 

Throughout  we  will  refer  to  the  generalized  shallow  water  equations  as 
with  “curvature”;  when  curvature  is  set  to  zero  it  is  termed  “no  curvature”,  and 
imposing  the  St.  Venant  equations  as  “  St.  Venant.”  The  test  input  is  shown 
in  Table  5.8.  The  model  has  the  same  geometry  and  node  layout  as  reported 
in  Section  3.5.  The  flow  rate  is  the  highest  tested  by  Sivakumaran.  Figure 
5.11  shows  a  comparison  of  water  surface  results  for  each  equation  set  and  the 
flume  data.  The  “curvature”  case  matches  the  flume  water  surface  quite  well. 
As  in  Chapter  3,  the  “no  curvature”  case  is  close  to  the  generalized  shallow 
water  equations  and  the  flume  data  in  the  supercritical  region  downstream  of 
the  crest,  but  shows  the  bulge  due  to  the  kinematic  inconsistency  we  have 
previously  discussed. 

The  water  surface  solution  for  the  St.  Venant  equations  starts  high 
and  drops  too  quickly.  This  difference  is  not  due  to  the  hydrostatic  assumption 
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since  the  “no  curvature”  case  is  also  hydrostatic.  Upstream,  the  St.  Venani 
equations  predict  a  higher  water  surface  to  maintain  this  flow  rate  but  without 
the  upstream  bulge  obtained  by  simply  setting  curvature  to  zero  in  the  gen¬ 
eralized  shallow  water  equations.  The  “no  curvature”  Ccise  bulge  in  a  concave 
region  is  due  to  an  overestimate  of  the  differential  volume.  This  in  turn  requires 
a  greater  pressure  gradient  to  accelerate  the  flow.  In  the  St.  Venant  equations, 
the  depth  is  measured  vertically  and  this  problem  is  relieved.  The  St.  Venant 
mild-slope  assumption  means  that  depth  is  measured  vertically  and  we  are  es¬ 
sentially  solving  the  horizontal  equations  with  a  modified  gravitation«il  force  to 
account  for  bed  slope. 

To  examine  this  mild-slope  cissumption  and  its  effects  upon  the  so¬ 
lution,  we  utilize  a  simple  comparison  of  the  St.  Venant  equations  (with  the 
mild-slope  assumption)  and  a  steep-siope  formulation.  The  steady-state,  one¬ 
dimensional,  steep-slope  shallow  water  equations  may  be  written: 


(5.4) 

(5.5) 


and  the  corresponding  St.  Venant  mild-slope  formulation  is: 


dx 


dx 


0 

0 


(5.6) 

(5.7) 


where,  Q  =  vd  and  P  =  uh.  All  other  geometric  terms  are  defined  in  Figure 
5.12.  Assuming  ^  to  be  a  constant,  the  steep-slope  formulation  can  be  written 
in  terms  of  the  variables  of  the  mild-slope  formulation  and,  after  solving  for 
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Figure  5.12:  Geometry  for  mild-slope  evaluati'^n. 
water  surface  slope  |^, 

dy  _  f  P^cos^g  1  db 

dx  \ghl~  P^cos'^dj  dx 

We  wish  to  compare  the  water  surface  slope  for  this  steep-slope  formulation  to 
that  of  the  St.  Venant  equations  near  h  —  hg.  The  mild-slope  result  is  obtained 
by  setting  ^  =  0.  There  is  a  slight  ambiguity  here  in  that,  =  tan  0,  but.  — 
is  actually  measured  and  is  not  among  the  terms  that  are  dropped  when  0  is 
assumed  to  be  small. 

First  we  consider  subcriticai  flow.  For  this  case  the  denominator  is 
positive  and  we  rewrite  (5.8) 

dy  _  f  P^cos^e  1^6 

dx  llghj- P^cos^e\l  dx  '  ' 


A  comparison  between  the  two  forms  can  be  made  by  setting  ^  =  0  in  the 
above  equation  for  the  mild-slope  form.  Both  have  water  surface  slopes  that 
are  in  the  opposite  direction  of  the  bed  slope.  It  is  also  apparent  that  the 
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Figure  5.13:  Subcritical  flow:  water  surface  slope  comparison. 


magnitude  of  is  larger  for  the  mild-slope  form.  Figure  5.13  shows  examples 
with  adverse  and  favorable  bed  slopes,  and  illustrates  the  relationship  between 
the  water  surface  slope  predicted  by  the  two  formulations. 


For  supercritical  flow  the  water  surface  slope  is  found  to  be 


dx 


■{ 


P^cosH 


Wo 


P^cos* 


,1 


dx 


(5.10) 


Here  the  mild-slope  assumption  again  exaggerates  the  water  surface  slope  but 
in  the  same  direction  as  the  bed  slope.  The  same  two  examples  under  super¬ 
critical  conditions  are  shown  in  Figure  5.14.  The  mild-slope  assumption  used 
in  the  St.  Venant  equations  tends  to  exaggerate  water  surface  slope  if  the  bed 
slope  is  steep.  The  amplification  is  roughly  related  to  1/  times  the  bed 

slope.  Hence  the  mild-slope  assumption  is  responsible  for  the  differences  on 
the  downstream  spillway  face  in  this  example,  not  curvature  effects. 


A  comparison  of  bed  pressures  is  shown  in  Figure  5.15.  The  gener¬ 
alized  shallow  water  equations  (curvature  case)  do  a  very  good  job  of  mod¬ 
eling  these  pressures  as  expected,  although  their  result  is  slightly  high  near 
X  =  —0.5m.  This  is  where  the  shallowness  parameter,  (curvature  multiplied  by 
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Figure  5.16:  Spillway  form:  shallowness  parameter. 


depth),  is  the  greatest.  (It  reaches  a  value  of  about  0.42,  see  Figure  5.16.) 
The  limit  given  by  Sivakumaran  on  the  applicability  of  these  equations  is 
—2<Kh<  0.54.  While  0.42  is  certzunly  within  these  limits,  we  must  note 
that  the  perturbation  analysis  implies  that  the  water  surface  “looks"  like  the 
bed;  i.e.,  the  flow  is  essentially  parallel  to  the  bed.  The  bed-normal  acceler¬ 
ations  become  more  important  as  this  is  relaxed.  These  effects  can  be  acco¬ 
modated  better  by  going  to  higher-order  terms  in  the  perturbation  analysis  as 
long  as  the  shallowness  parameter  is  not  too  large.  This  is  a  serious  limitation 
to  applicability  of  the  model  to  a  large  class  of  spillways  which  typically  have 
a  steep  upstream  face.  Downstream,  the  slight  discrepancy  near  the  toe  of  the 
spillway  has  been  attributed  by  Sivakumaran  to  the  irrotational  flow  assump¬ 
tion  (over  depth).  This  is  based  upon  the  research  of  Henderson  and  Tierney 
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[28].  Overall,  the  comparison  is  quite  remarkable.  On  the  other  hand,  the  St. 
Venant  and  the  “no  curvature”  cases  are  poor.  This  is  not  suprising  since  both 
are  hydrostatic  and  the  St.  Venant  model  measures  pressure  based  on  depth 
measured  vertically. 

These  results  indicate  that  the  longitudinal  profile  of  water  surface 
and  bed  pressure  are  modeled  fairly  well  with  the  generalized  shallow  water 
equations.  The  new  model  is  limited  in  handling  some  common  spillways  in 
which  the  upstream  face  is  steep.  Downstream  of  the  crest,  where  the  flow'  is 
supercritical,  the  water  surface  more  closely  parallels  the  bed  and  thus  the  per¬ 
turbation  analysis  is  more  accurate.  The  “no  curvature”  case  of  the  generalized 
shallow  water  equations  also  compares  well  with  the  flume  water  surface.  This 
implies  that  the  hydrostatic  pressure  assumption  is  reasonable  in  the  down¬ 
stream  supercritical  reach.  The  overall  poor  comparison  with  the  St.  Venant 
model  can  be  attributed  to  the  mild-slope  assumption.  The  poor  spillway  ca¬ 
pacity  prediction  is  a  result  of  the  hydrostatic  assumption  in  the  steep-slope 
and  the  St.  Venant  models. 

5.2.3  Outletworks  Flume. 

The  final  test  is  the  most  general  we  shall  undertake.  We  compare 
the  previous  sets  of  model  equations  to  an  outletworks  physical  model  at  the 
Waterways  Experiment  Station,  WES.  Here  our  primary  concern  is  the  lateral 
behavior  predicted  by  the  mathematical  models  for  more  realistic  structures. 
A  diagram  of  the  flume,  as  tested,  is  shown  in  Figure  5.17.  The  plan  view 
geometry  for  the  flume  and  the  model  are  identical  but  the  numerical  approx- 
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Figure  5.17:  Outletworks  flume:  flume  and  numerical  model  geometry. 
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imation  to  bed  elevation  has  smooth  transitions  to  provide  a  continous  slope. 
The  dimensions  of  the  flume  were  chosen  such  that  typical  operating  condi¬ 
tions  place  the  flow  regime  within  the  fully  turbulent  region.  For  our  test,  the 
Reynolds  number  is  about  ft  =  ~  ss  200000,  which  indicates  that  scale-effects 
are  small.  For  even  relatively  smooth  bed  and  wall  conditions  the  flume  will 
be  in  the  fully  turbulent  regime  (see  (48,  13]).  Therefore,  w-  might  reasonably 
expect  the  flume  results  to  exhibit  behavior  similar  to  large-scale  structures  of 
this  type. 

The  bed  elevation  of  the  flume  is  purposely  laid  out  longitudinally 
as  a  parabolic  function  to  minimize  flow  separation.  The  actual  dimensions  in 
Figure  5.17  are  “as  measured”  and  vary  slightly  from  the  construction  spec¬ 
ifications,  primarily  as  a  result  of  settling.  The  transition  from  the  circular 
conduit  to  the  apron  is  fairly  complicated.  The  walls  have  a  slight  flair  with  a. 
radius  of  2.75ft  and  the  bed  has  fillets  to  move  smoothly  from  the  round  con¬ 
duit  to  the  rectangular  section  of  the  apron.  The  apron  wall  has  an  outward 
slope  of  approximately  g.  A  level  stilling  basin  which  contains  baffles  and  an 
end  sill  is  located  4.015ft  from  the  conduit  exit.  These  stilling  basin  features 
are  not  numerically  modeled  here;  since  these  effects  far  downstream  are  of 
lesser  importance  to  the  spillway  flow.  The  actual  bed  elevation  numerically 
modeled  (as  shown  in  Figure  5.17)  contains  no  discontinuities  in  slope  so  that 
the  curvature  exists  throughout.  The  bed  is  flat  in  the  lateral  direction,  i.e., 
all  slope  and  curvature  variations  are  in  the  longitudinal  direction. 

The  flume  water  surface  is  recorded  every  0.05ft  laterally  and  every 
0.10ft  longitudinally.  An  additional  reading  is  made  near  the  sidewalls  (typi¬ 
cally  within  0.02ft).  The  measurements  are  recorded  with  a  point  gage  capable 
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of  measuring  to  a  precision  of  0.001ft  (it  is  graduated  in  0.01  ft  units  and  with 
a  vernier  one  can  detect  0.001  changes).  In  reality  the  water  surface  is  quite 
rough  and  time  varying,  so  that  readings  are  approximately  accurate  to  about 
0.02ft.  The  flow  rate  is  measured  using  a  12  by  6  inch  venturi  meter  with  a 
differential  manometer.  The  inflow  pressure  is  maintained  by  a  constant  head 
tank.  The  manometer  calibration  is  checked  using  direct  measurements  of  vol¬ 
ume  change  in  a  specific  time  period.  The  discharge  is  estimated  [57]  to  be 
accurate  within  3  %. 

Tht  flume  water  surface  results  are  shown  in  Figure  5.18  as  a  three- 
dimensional  perspective  for  a  flow  rate  of  2.56  cfs  (cubic  feet  per  second).  The 
view  shows  the  lateral  hump  in  the  water  surface  as  the  flow  exits  the  outlet 
conduit.  This  eventually  spreads  to  a  relatively  uniform  depth  as  intended  in 
the  design.  If,  however,  flow  becomes  focussed  the  result  will  be  circulation  in 
the  stilling  basin  with  possible  damage  there.  The  nonuniform  velocity  field 
would  also  cause  more  downstream  erosion.  The  drag  of  the  apron  sidewalls 
causes  the  sharp  rise  in  water  surface.  The  buildup  in  water  depth  along  the 
sides  of  the  stilling  basin  is  a  result  of  the  change  in  side  slopes. 

The  numerical  model  is  constructed  of  408  nodes  and  368  elements 
(see  Figure  5.19).  It  extends  longitudinally  from  x  =  0.75ft  to  i  =  6.0ft. 
The  upstream  boundary  is  therefore  below  the  complex  geometry  of  the  fil¬ 
lets  near  the  conduit  exit  and  the  vertical  velocity  profile  should  be  closer  to 
an  open-channel  flow  distribution.  The  flow  is  supercritical  throughout  the 
study  reach,  and  at  the  upstream  boundary  the  velocity  and  depths  are  speci¬ 
fied.  The  depths  are  the  flume  results  interpolated  from  the  0.70ft  and  0.80ft 
ranges.  The  x  component  velocity  distribution  is  assumed  to  be  uniform  across 


Figure  5.18:  Outletworks  flume:  water  surface  results,  perspective  view. 
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Figure  5.19:  Outletworks  model:  computational  mesh. 


the  upstream  boundary.  The  lateral  velocity  component  at  the  inlet  is  speci¬ 
fied  so  that  the  flow  is  tangent  to  the  sidewadls  and  zero  in  the  channel  center 
with  linear  variation  between  these  stations.  These  are  the  simplest  reasonable 
boundary  conditions  and  should  be  suitable  for  a  comparison  of  model  equa¬ 
tions.  Furthermore,  this  is  as  much  information  as  a  modeler  typically  hea  to 
make  a  site-specific  study.  Other  model  conditions  are  shown  in  Table  5.9.  The 


Table  5.9:  Outletworks  model. 


Input  Parameters 

9 

n 

0.012 

Vt 

0.05ft^/sec 

At 

O.OSsec 

Simulation  Time 

I5sec 

Iterations 

3 

t 

0.03 
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Figure  5.20:  Outletworks  flume:  depth  contours. 

turbulent  viscosity  corresponds  to  a  value  of  Cb  of  about  O.I. 

Figures  5.20  -  5.23  show  the  depth  contours  for  the  flun;’.  the  general¬ 
ized  shallow  water  equations,  the  standard  steep-slope  shallow  water  equations, 
and  the  St.  Venant  equations.  The  flume  and  St.  Venant  measurements  are 
made  vertically  as  water  surface  elevations  which  we  have  converted  to  depths 
normal  to  the  bed.  The  flume  data  in  the  upper  region  reveal  a  cop*’nuation  of 
the  conduit  shape,  i.e.,  a  hump  in  the  water  surface  with  the  centerline  being 
a  relative  maximum.  This  form  ends  near  i  =  1 .75ft.  From  there  on  the  cen¬ 
terline  is  a  relative  minimum  or  the  surface  is  nearly  flat.  This  is  the  desired 
pattern  for  a  good  hydraulic  design.  At  the  channel  edges  the  wall  drag  causes 
a  substantial  buildup  of  water  depth. 

A  most  important  result  of  the  generalized  shallow  water  equations 
is  the  more  gradual  lateral  spread  of  the  flow  tube.  As  an  illustration,  within 
the  portion  of  the  flume  modeled  here,  the  flume  shows  the  centerline  to  be 
a  relative  maximum  for  about  1.0ft.  The  model  including  curvature  gave  a 
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Figure  5.21;  Outletworks  model:  generalized  shallow  water  equations,  depth 
contours. 


Figure  5.22;  Outletworks  model:  standard  steep-slope  shallow  water  equations, 
depth  contours. 
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Figure  5.23:  Outletworks  model:  Si.  Venant  equations,  depth  contours. 


similar  length.  The  steep-slope  shallow  water  model  (no  curvature)  and  the 
St,  Venant  equations  are  nearly  identical  throughout  the  domain.  The  length 
for  which  the  centerline  is  a  maximum  as  predicted  by  these  models  is  only 
about  0.7ft.  Hence,  the  nonhydrostatic  component  has  a  definite  effect  upon 
the  wave  speed  here.  This  is  a  region  in  which  the  flow  is  entirely  supercritical 
and  for  which  our  one-dimensional  test  showed  little  improvement  by  includ¬ 
ing  the  nonhydrostatic  pressure  component  in  the  generalized  shallow  water 
equations.  However,  now  the  lateral  variation  is  strongly  dependent  upon  the 
nonhydrostatic  contribution.  The  predicted  wave  celerity  for  a  similar  case  in 
which  the  curvature  and  dominant  flow  direction  are  along  the  x  axis  can  be 
shown  (see  [20])  to  be 


log(l  -  Kyh) 

(1  -  Kih)Ki 


On  the  other  hand  the  celerity  without  curvature  included  is  simply 


I 

^S3we  ^  (93^)^ 


(5.12) 
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OUTLETWORKS  FLUME 
LATERAL  DEPTH  VARIATION 


DISTANCE .  f=T 

Figure  5.24:  Outletworks  flume;  lateral  depth  profiles. 

A  convex  curvature  (kj  negative)  produces  a  lower  celerity  in  the  generalized 
equations  which  is  not  captured  in  the  standard  shallow  water  equations.  The 
true  solution,  which  includes  nonhydrostatic  pressure  contributions  as  a  result 
of  short  wavelengths  in  addition  to  bed  curvature  effects,  has  a  wave  celerity 
that  is  even  smaller  than  the  generalized  equations  predict.  In  a  concave  region 
(positive  curvature)  the  wave  celerity  predicted  by  the  generalized  shallow  water 
equations  will  be  greater  than  from  the  standard  shallow  water  equations. 

Lateral  profiles  of  depth  are  shown  in  Figures  5.24-5.26.  The  fig¬ 
ures  include  flume,  generalized  shallow  water,  and  standard  steep-slope  shallow 
water  results.  The  St.  Venant  results  are  nearly  identical  to  the  standard  steep- 
slope  results.  The  lateral  profiles  are  at  distances  downstream  of  the  conduit 
of  1.0,  1.6,  2.8,  and  5.0ft.  These  results  provide  an  additional  view  of  the  infor- 
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DEPTH 


102 


OUTLETWORKS  MODEL  (NO  CURVATURE) 
UTERM.  DEPTH  VARIATION 


DISTANCE .  FT 


Figure  5.26;  Outletworks  model:  standard  steep-slope  shallow  water  equations, 
lateral  depth  profiles. 
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mation  contained  in  the  contour  plots.  Between  1.0ft  and  1.6ft  the  flume  and 
model  including  curvature  effects  show  a  gradual  collapse  of  the  wave.  With 
no  curvature  effects  the  wave  has  propagated  out  to  the  edges  leaving  a  trough 
in  the  center.  Downstream,  the  distinction  between  the  results  from  the  two 
model  equation  sets  is  very  small.  The  model  in  both  cases  only  captures  the 
behavior  of  the  flume  in  a  qualitative  sense.  One  reason  is  that  we  use  actual 
flume  measurements  for  the  water  surface  upstream  boundary  condition,  but 
rely  upon  a  uniform  velocity  profile.  A  second  reason  is  that  the  wavelength  of 
the  hump  in  the  chajinel  center  is  only  about  1ft  while  the  the  depth  is  nearly 
0.4ft.  This  collapse  rate  is  really  a  short-wave  phenomenon  not  capable  of  be¬ 
ing  accurately  represented  by  a  long-wave  model.  However,  the  effects  of  the 
bed  curvature  are  captured  and  do  improve  the  comparison  in  the  generalized 
shallow  water  model. 

The  lateral  velocity  profile  comparisons  between  the  generalized  shal¬ 
low  water  equations  and  the  standard  steep-slope  shallow  water  equations  are 
shown  in  Figures  5.27-5.30.  The  distinctive  diff'erences  in  these  figures  occur 
at  X  =  1.0ft  and  i  =  2.8ft.  At  x  =  l.Oft  the  pressure  is  lower  for  the  case 
with  curvature,  and  since  velocity  is  dependent  upon  the  pressure  gradient,  the 
lateral  velocity  is  lower.  At  i  =  2.8ft,  the  model  with  curvature  produces  gen¬ 
erally  higher  lateral  velocities.  The  outward  propagation  of  the  lateral  wave  is 
quicker  in  the  “no  curvature”  case  and  the  peak  has  propagated  out  to  the  wall, 
whereas,  for  the  case  including  curvature,  the  outer  depth  has  not  reached  its 
peak.  The  result  is  a  reduced  adverse  gradient  at  this  point  in  the  generalized 
equations,  and  thus  slightly  larger  outward  velocity. 

The  mild-slope  assumption  made  little  difference  for  this  case;  i.e.,  the 
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COMPARISON  OF  LATERAL  VELOCITY 


Xa  1.0FT 


DISTANCE ,  FT 

Figure  5.27:  Outletworks  flume:  lateral  velocity  profiles,  x  =  1.0ft. 

St.  Venant  equations  yielded  results  quite  close  to  those  of  the  stauidard  steep- 
slope  equations.  On  the  other  hand  the  bed  curvature  effects  are  important  and 
the  generalized  shallow  water  equations  did  give  better  results  for  the  water 
surface  distribution.  However,  the  primary  concerns  around  an  outletworks  are 
related  to  the  spread  of  the  flow  tube  which  is  a  short-wave  phenomena  and 
none  of  the  models  performed  well  in  this  regard. 


5.3  Discussion 

The  most  striking  contrast  between  results  of  the  generalized  shallow 
water  equations  and  the  other  formulations  is  in  the  predictions  of  spillway 
capacity.  This  is  indeed  a  result  of  the  hydrostatic  pressure  distribution  used 
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in  the  standard  equations.  In  reality  the  negative  curvature  near  the  crest  of  the 
spillway  reduces  the  pressure  and  so  allows  more  flow  for  a  given  water  level. 
We  shall  develop  an  understanding  of  this  by  considering  a  linearization  of 
the  generalized  shallow  water  equations  that  will  allow  us  to  estimate  spillway 
capacity  for  the  standard  methods  and  the  change  expected  by  implementation 
of  the  generalized  shallow  water  equations. 

Consider  the  discharge  per  unit  width  (for  a  one-dimensional  flume) 
as  defined  in  the  generalized  shallow  water  equations  to  be  given  by 

qu;=  u(s3)dsj  (5.13) 

Jo 

or, 

=  (5-14) 

where,  Uh  is  the  velocity  at  the  surface  and  as  before  f  {h)  =  1  —  k/i.  The 
specific  energy  is  given  by 

E.  =  $  +  h  (5.15) 

This  is  the  energy  at  the  surface  with  the  datum  at  the  spillway  crest.  We  wish 
to  compare  the  spillway  capacity  as  determined  by  the  specific  energy  with 
curvature  and  for  no  curvature.  The  “no  curvature’’  ceise  corresponds  to  the 
St.  Venant  model  and  the  standard  steep-slope  shallow  water  equations.  For  a 
given  discharge  the  specific  energy  is  a  minimum  at  the  spillway  crest.  We  can 
see  this  by  considering  an  idealized  flow  in  a  subcritical  channel.  Across  the 
channel  is  a  “bump”  that  is  originally  of  very  small  height.  Now  as  we  gradually 
raise  the  height  of  the  bump  a  choke  point  is  eventually  reached  where  any 
further  increase  in  the  “bump”  height  will  cause  the  upstream  water  surface  to 


rise.  Additional  energy  is  required  to  pass  this  flow  for  any  additional  increase 
in  height.  Also,  downstream  of  the  bump  the  flow  is  no  longer  subcritical.  So 
at  the  crest  of  the  bump  the  flow  is  termed  critical  and  the  specific  energy  is  a 
minimum. 

Now  using  equations  5.14  and  5.15  for  x  —  and  |  x  l<<  where 
ho  is  a  nominal  value  of  h,  an  approximation  for  specific  energy  is 

£a«^(l  +  x)  +  /^  (5.16) 

Minimizing,  we  find  ho,  the  critical  depth 

—  (1  +  X)  (5.17) 

9 

or  in  terms  of  specific  energy 

£.«lf4(‘+x)r  (5-18) 

2  L  9 

This  demonstrates  that  if  x  <  0,  as  it  would  be  near  the  crest  the  specific 
energy  is  less  than  for  the  no  curvature  case,  x  =  0.  Therefore,  we  would 
expect  the  St.  Venant  model  and  other  models  not  including  bed  curvature  to 
underestimate  spillway  capacity. 

If  we  solve  for  the  discharge  in  terms  of  the  required  specific  energy 
we  obtain  the  spillway  capacity.  The  ratio  qr  of  discharge  without  curvature 
terms  to  that  with  curvature  included  is  approximated  by 

9r«5(l  +  x)’  (5.19) 


as  demonstrated  in  Figure  5.31. 
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Hence  we  see  that  near  the  crest  the  curvature  terms  are  very  impor¬ 
tant.  Models  not  including  this  nonhydrostatic  effect  will  overestimate  critical 
depth  and  underestimate  spillway  capacity.  As  shown  in  Figure  5.31  this  can 
be  significant.  The  two-dimensional  model  can  treat  this  along  with  the  side 
constrictions  to  capture  the  spillway  capacity  more  accurately. 
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Figure  5.29:  Outletworks  flume:  laterd  velocity  profiles,  x  =  2.8f1 
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Figure  5,30:  Outletworks  flume:  lateral  velocity  profiles,  i  =  5. 
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Chapter  6 


Conclusion 


In  this  investigation  we  have  developed  a  generalized  shallow  water 
equation  set  to  reproduce  free-surface  flow  over  curved  beds  and  applied  the 
model  to  spillway  flows.  These  generalized  shallow  water  equations  include 
bed  curvature  effects,  are  nonhydrostatic,  and  are  two-dimensional.  The  equa¬ 
tions  are  derived  via  a  perturbation  expansion  in  a  shallowness  parameter.  A 
significant  assumption  made  in  this  derivation  is  that  the  flow  is  irrotational 
about  axes  parallel  to  the  plane  of  the  bed.  However,  there  is  no  limitation  on 
vorticity  about  axes  normal  to  the  bed. 

A  one-dimensional  model  was  constructed  to  assess  numerical  prob¬ 
lems  and  determine  additional  needs  for  the  final  two-dimensional  investigation. 
This  finite  element  model  w2is  straightforward  and  employed  time- lagging  of 
nonlinear  terms.  We  considered  available  flume  data  for  a  simple  curved  bed 
profile.  In  comparing  the  one-dimensional  equations  with  and  without  curva¬ 
ture  the  principal  differences  occured  upstream  of  the  spillway  crest.  This  is 
important  since  this  is  the  spillway  capacity. 

The  numerical  handling  of  the  convection  terms  used  a  method  of 
artificial  viscosity  weighted  by  the  convection  velocity.  At  first  a  Taylor  Weak 
Form  approach  was  attempted  but  was  shown  to  have  some  difficulty  near 
points  of  transition  from  supercritical  to  subcritical  regimes  and  vice  versa. 
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The  artificial  viscosity  avoided  these  problems  but  also  has  some  detractions 
and  was  replaced  in  the  final  two-dimensional  model.  (One  obvious  problem 
with  the  viscosity  approach  is  adding  diffusion  terms  to  the  mass  conservation 
equation.) 

The  final  development  involved  a  two-dimensional  model  in  which  the 
nonlinear  terms  were  addressed  via  a  Newton  Raphson  technique,  iterating  at 
each  time  step.  This  allows  larger  time  steps.  We  introduced  a  Petrov  Galerkin 
scheme  based  on  the  SUPG  idea  and  A-schemes  in  which  the  test  function  is 
weighted  using  the  eigenvalues  of  the  convection  matrices. 

A  two  stage  testing  program  was  conducted.  The  first  stage  vali¬ 
dated  the  discrete  model  against  analytic  benchmark  cases  and  the  second 
stage  compared  the  model  results  to  flume  data.  The  first  stage  was  accom¬ 
plished  by  several  one-dimensional  tests  performed  along  the  si  and  then  along 
S2  directions.  These  tests  were  for  both  constant  and  variable  curvature.  The 
Bernoulli  equation  served  as  the  standard  for  comparison.  This  validated  the 
discrete  model  for  the  derived  equations. 

Tests  of  the  suitability  of  the  equations  were  then  undertaken  by 
comparison  to  actual  flume  data.  The  first  study  was  a  comparison  to  a  su¬ 
percritical  transition  in  which  there  was  no  bed  curvature.  This  was  primarily 
a  test  of  the  numerical  scheme  and  boundary  condition  implementation.  The 
model  performed  well  within  the  limitations  of  the  shallow  water  equations. 

The  next  set  of  tests  was  made  for  a  curved  bottom  flume  with  straight 
walls.  Results  from  the  generalized  shallow  water  equations  compared  remark¬ 
ably  well  to  the  flume  water  surface  data.  The  spillway  capacity  was  accurately 
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captured  and  a  comparison  of  bed  pressure  was  quite  good.  The  St.  Venant 
and  the  standard  steep-slope  shallow  water  equations  missed  the  spillway  ca¬ 
pacity,  both  prediciting  a  greater  water  surface  elevation  to  attain  the  steady 
discharge.  (The  hydrostatic  pressure  aissumption  nnisses  the  reduction  in  pres¬ 
sure  caused  by  the  uplift  of  the  flow  over  the  curved  crest.  This  is  responsible 
for  the  error  in  spillway  capacity.) 

Along  the  downstream  face  the  standard  steep-slope  formulation  com¬ 
pared  fairly  well  with  the  generalized  shallow  water  equations  and  flume  data. 
The  St.  Venant  equations  predicted  too  steep  a  water  surface  slope.  This  was 
shown  to  be  a  result  of  the  mild-slope  assumption,  and  not  curvature  related. 
The  principal  improvement  of  the  new  equations  over  the  ‘‘no  curvature  case" 
is  in  the  neighborhood  of  the  spillway  crest. 

We  collected  experimental  data  on  a  flume  at  the  Waterways  Exper¬ 
iment  Station  designed  to  study  outlet  works.  Here  the  emphzisis  was  upon 
lateral  variations  in  water  surface.  The  improvement  offered  by  the  new  equa¬ 
tions  was  subtle  showing  slower  (and  more  realistic)  spread  of  the  flow  conduit 
over  negative  curvature  regions.  Overall,  however,  the  representation  of  what  is 
essentially  a  short  wave  phenomena  was  weak.  However,  the  results  did  provide 
a  qualitative  picture  of  the  flow  that  is  beneficial  for  designers. 

A  list  of  the  conclusions  then  is  as  follows: 

1 .  These  equations  are  a  significant  improvement  over  standard  steep-slope 
shallow  water  equations  in  the  vicinity  of  the  spillway  crest.  Within  the 
stated  restrictions  on  the  size  of  curvature  and  depth,  these  equations 
can  provide  a  good  representation  of  spillway  capacity. 
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2.  These  equations  do  show  an  improvement  in  the  lateral  depth  distribu¬ 
tion  though  they  are  shallow  water  equations  and  suffer  in  predicting 
phenomena  associated  with  short  wavelengths. 

3.  Along  the  downstream  face  of  the  spillway  or  outlet  channel  the  flow  is 
supercritical  and  the  difference  in  water  surface  attributable  to  curvature 
is  small,  generally  on  the  order  of  precision  used  in  physical  model  flumes. 

4.  The  St.  Venant  equations  proved  to  be  much  worse  than  the  generalized 
or  the  standard  steep-slope  shallow  water  equations  along  the  downstream 
face.  The  difference,  however,  is  attributable  to  the  mild-slope  assump¬ 
tion  and  not  to  the  hydrostatic  pressure  «issumption  or  other  curvature 
manifestations. 


Recommendations  for  Future  Work 

The  equations  are  most  important  in  the  region  around  the  spillway 
crest  and  for  evaluation  of  contraction  effects  near  the  ends  of  the  spillway. 
This  is  where  the  nonhydrostatic  pressure  distribution  is  most  critical.  Future 
research  should  concentrate  on  this  issue.  Furthermore: 

1.  Many  spillways  have  steep  faces  and  abrupt  curvature  that  violate  the 
limits  imposed  by  the  shallowness  condition  in  this  perturbation  expan¬ 
sion.  This  significantly  reduces  the  pressure  at  the  crest  and  so  increases 
spillway  capacity.  To  be  a  practical  tool,  the  model  needs  to  be  extended 
to  this  more  general  case,  perhaps  by  the  calculation  of  a  flow  line  repre¬ 
senting  the  bed. 
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2.  Higher-order  terms  in  the  expansion  for  the  steep-slope  shallow  water 
equations  can  also  be  constructed  for  this  region.  In  particular,  do  these 
equations  represent  this  region  as  well  as  the  generalized  shallow  water 
equations?  Of  course,  these  Boussinesq  equations  include  third-order 
spatial  derivatives  and  so  will  require  additional  computational  effort  and 
some  improved  convection  scheme. 

3.  The  unsteady  behavior  of  the  equations  could  be  quite  interesting.  In 
fact,  even  for  steady  boundary  conditions  the  highly  nonlinear  flow  near 
and  along  the  spillway  face  is  unsteady. 


Appendix 


Integration  of  Coefficients 

The  integration  of  the  coefficients  used  in  the  generalized  shallow 
water  equations  was  accomplished  using  the  software  Mathematica  [43]  and  are 
as  follows: 
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